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ABSTRACT 

Context. With the high number of extrasolar planets discovered by now, it becomes possible to use the properties of this planetary 
population to constrain theoretical formation models in a statistical sense. This paper is the first in a series in which we carry out a 
large number of planet population synthesis calculations within the framework of the core accretion scenario. We begin the series 
with a paper mainly dedicated to the presentation of our approach, but also the discussion of a representative synthetic planetary 
population of solar like stars. In the second paper we statistically compare the subset of detectable planets to the actual extrasolar 
planets. In subsequent papers, we shall extend the range of stellar masses and the properties of protoplanetary disks. 
Aims. The last decade has seen a large observational progress in characterizing both protoplanetary disks, and extrasolar planets. 
Concurrently, progress was made in developing complex theoretical formation models. The combination of these three developments 
allows a new kind of studies: The synthesis of a population of planets from a model, which is compared with the actual population. 
Our aim is to get a general overview of the population, to check if we quantitatively reproduce the most important observed properties 
and correlations, and to finally make predictions about the planets that are not yet observable. 

Methods. Based as tightly as possible on observational data, we have derived probability distributions for the most important initial 
conditions for the planetary formation process. We then draw sets of initial conditions from these distributions and obtain the corre- 
sponding synthetic planets with our formation model. By repeating this step many times, we synthesize the populations. 
Results. Although the main purpose of this paper is the description of our methods, we present some key results: We find that the 
variation of the initial conditions in the limits occurring in nature leads to the formation of planets of large diversity. This formation 
process is best visualized in planetary formation tracks in the mass-semimajor axis diagram, where different phases of concurrent 
growth and migration can be identified. These phases lead to the emergence of sub-populations of planets distinguishable in a mass- 
semimajor axis diagram. The most important ones are the "failed cores", a vast group of core-dominated low mass planets, the 
"horizontal branch", a sub-population of Neptune mass planets extending out to 6 AU, and the "main clump", a concentration of giant 
gaseous giants planets at around 0.3-2 AU. 

Key words. Stars: planetary systems - Stars: planetary systems: formation - Stars: planetary systems: protoplanetary disks - Planets 
and satellites: formation - Solar system: formation - Methods: numerical 



1. Introduction 

As of spring 2009, more than 300 extrasolar planets have been 
discovered (J. Schneider's Extrasolar Planet Encyclopedia at 
http://exoplanet.eu). The richness and diversity of the charac- 
teristics of these exoplanets like their mass or semimajor axis 
is impressive, and was not necessarily expected from the single 
example - our own solar system - that was available to study 
before t he dis covery of 51 Peg b (HD217014b) by Mayor & 
Queloz JT9951 

Since then, the observational field of extrasolar planet search 
has seen a rapid evolution leading to numerous additional dis- 
coveries of planets orbiting other stars. These discoveries have 
also triggered numerous theoretical studies about the formation 
and evolution of these planets. Key physical processes in planet 
formation and evolution could be identified whose importance 
was not fully realized in previous works based on the solar sys- 
tem alone. 
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Some of these discovered planets, and multiple planetary 
systems, are sufficiently interesting by themselves to warrant in- 
dividual theoretical studies. Examples are the extrasolar plane- 
tary system with three Neptune-mass planets around HD 69830 
(Lovis et al. 120061 Alibert et al. 12006b . or the transiting Neptune 
mass planet GJ436 b (Butler et al. 120041 Gillon et al. [20071 
Figueira et al. 2008 ). Of course, the giant planets in our own so- 
lar system provide a much larger and detailed set of constraints 
than any known extrasolar planet. Therefore, each formation 
model applied to discuss extrasolar planet formation should also 
be put to the test to reproduce the characteristics of our own gi- 
ant planets (Pollack et al. 1 19961 Alibe rt et al. l2005bl Hubickyj et 
al. 120051 Benvenuto & Branini 120051 

The modeling of the formation of such single systems while 
a necessary condition to validate formation models is not sat- 
isfactory by itself. Indeed, the number of model parameters is 
generally large while the number of constraints deriving from 
a single system is small, and not strong enough to completely 
constrain any formation model. 

Thanks to the rapid growth of the number of known extra- 
solar planets, the situation has however dramatically changed: 



2 



C. Mordasini et al.: Extrasolar planet population synthesis I 



Instead of having only a single object or a single system to study, 
we now begin to be able to describe an entire population of ex- 
trasolar planets orbiting FGK stars in the solar neighborhood. 
While this population is still smaller than one would ideally like, 
it nevertheless already allows to extract statistically a wealth of 
information (e.g. Udry & Santos 2007; Cumming et al. 2008) 
to constrain formation models that exceeds by far what one ex- 
trasolar planet can do. This is especially true since most of the 
extrasolar planets have been discovered by radial velocity mea- 
surements so that only a few orbital elements and a minimum 
mass are known for one individual object. For the growing num- 
ber of transiting planets more physical properties can be derived 
and compared with internal structure models (Baraffe et al. 2008 ; 
Figueira et al. 2008 ). Unfortunately, transiting planets known so 
far are all in close proximity to their host star. Hence it is some- 
times unclear to what extend their characteristics are still related 
to their formation or rather to subsequent evolution (e.g. evapo- 
ration). 

Parallel to the discovery of more and more end-products of 
the planetary formation process i.e. planets, large observational 
progress (e.g. Meyer et al. 12006 ) has also been made in charac- 
terizing the initial conditions for this process, i.e. the protoplan- 
etary disks. Thanks to these observations, we begin to be able to 
determine the probability of occurrence of any particular initial 
condition for planetary formation, like disk metallicity, mass or 
lifetime. 

With these two sets of observational data at hand, a new in- 
teresting class of theoretical planet formation studies has become 
possible, where a theoretical model serves as the link between 
these two groups of observations: The synthesis of populations 
of planets by Monte Carlo methods. In this approach the ob- 
served distributions of disk properties are used as varying initial 
conditions for the model. The final characteristics of the syn- 
thetic planets that form in the model can then be compared sta- 
tistically to those of the actual observed populations. This ad- 
dresses the question if the observed diversity of extrasolar plan- 
ets is simply the consequence of the diversity of disk properties. 

As we shall show, such studies have proven to be very fruit- 
ful, as they not only allow to reproduce observations but also 
show the links and correlations between the different initial con- 
ditions and the characteristics of the resulting planets. Thereby 
they provide great insights into the formation mechanism. Last 
but not least, such an approach by predicting the actually existing 
planet population as opposed to the actually detected one, allows 
to optimize future searches and instruments when coupled to a 
synthetic detection bias for a particular detection method. 

Compared to similar studies e.g. the pioneering work of Ida 
& Lin <2004all2004bll2005l a ndl2D^ . or the studies of Kornet& 
Wolf d2006l Robinson et al. d2006t and Thommes et al. d2008l . 
we have attached particular importance to three distinct areas: 
First, we use the detailed extended core accretion formation 
model of Alibert et al. (2005a) that has been successfully ap- 
plied to quantitatively explain the many observed constraints of 
the giant planets of our own solar system (Alibert et al. [2005b ). 
Second, we stick as tightly as possible to probability distribu- 
tions derived from observations, and third (cf. the companion 
paper Mordasini et al. (12008 >, hereafter paper II), we use quan- 
titative statistical methods to compare model outcomes and ob- 
servations and require that as many different observational con- 
straints as possible be satisfied at the same time. In this way, we 
can check which observed properties can be reproduce by our 
formation model. But also discrepancies between the synthetic 
and the actual population provide new insights allowing to im- 
prove the models. 



In this first aper, we present the methods we use to generate 
the synthetic population, in particular the formation model and 
the probability distributions for the Monte Carlo variables. We 
then show the resulting planetary formation tracks in the semi- 
major vs. mass plane. These tracks are of fundamental impor- 
tance to understand the characteristics of the resulting synthetic 
population, as they illustrate how and why planets reach their 
final position in the mass-distance diagram. This a- M distribu- 
tion at the end of the formation phase is characterized by a num- 
ber of structures (clumps, concentrations and depletions). Some 
particular regions in this diagram are identified and discussed. 
In the companion paper II, we use a synthetic detection bias for 
the radial velocity method to identify the subset of detectable 
synthetic planets. We then compare this sub-population with 
statistical, quantitative methods to the actual extrasolar planet 
population. In these first two papers, we assume a mass of the 
host star of 1 M Q , as most known extrasolar planets are found 
around solar type stars. In later papers in this series we shall 
study the influence of different host star masses as well as of 
different formation environments (i.e. disk properties). Some of 
these have already been observationally identified such as the 
well known correlation between the stellar metallicity and the 
detection probability of giant planets. 

The outline of the paper is as follows: In section ^2]we give 
an overview of our formation model, with a focus on the modi- 
fications and necessary simplification^] relative to Alibert et al. 
(2005a). In section ^3] the Monte Carlo approach is described, 
whereas we determine the probability distributions for the initial 
conditions in §|4] As results, section ^5] illustrates the numeri- 
cal population synthesis process with example formation tracks 
in the distance-mass plane and discusses the properties of the 
planet population. The conclusions are drawn in the last section, 



2. Giant planet formation model 

The link between the initial conditions i.e. the properties of the 
protoplanetary disk, and the final outcome i.e. the planets, can 
only be given by a theoretical formation model. This link is usu- 
ally very complicated involving many feed-back mechanisms 
and nonlinearities. For the simulations presented in this paper, 
we calculate in a consistent way the formation of a protoplanet, 
its migration, and the structure and evolution of the protoplane- 
tary disk. 



2. 1 . Disk structure and evolution 

The structure and evolution of the protoplanetary disk is mod- 
eled as a non-irradiated, 1 + 1D c-disk (Shakura & Sunyaev 
1973), following the method originally presented in Papaloizou 
& Terquem ( 1999). We thus solve the diffusion equation (effec- 
tive viscosity v) describing the evolution of the gas surface den- 
sity 2 as a function of time t and distance a to the star: 



dt 



3d_ 
a da 



,1/2. 



da 



(via 1 ' 2 ) 



+ i,,(fl) 



(1) 



1 The synthesis of a population of ~ 30 000 planets takes several days 
on a 50 CPU cluster. Most of the time is spent in solving the planetary 
envelope structure equations. 
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The photo-evaporation term £„, is given by (Veras & Armitage 
2003): 
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for a < R g 
otherwise 



(2) 



where R g is taken to be 5 AU, a max is the size of the disk, and the 
total mass loss M w due to photo-evaporation is an input parame- 
ter which together with the a parameter determines the lifetime 
of the disk. 

For simplicity, we adopt an initial profile of the gas disk 
surface density according to the phenomenological model of 
Hayashi JlWTT) . £(a,f = 0) = S (a/a )" 3/2 where Z is the 
surface density at our reference distance (ao=5.2 AU), and the 
computational disk extends from a m i n =0.1 AU to a max =30 AU. 
The initial total gas disk mass in the computational disk is then 

47rEoflo /2 (amax - a ]^)- For tne initial profile with S oc aT 312 the 
accretion rate decreases from the inner to the outer parts of the 
disk. As shown in e.g. Papaloizou & Terquem (1999), the in- 
ner parts of the disk evolve rapidly toward a state of constant 
accretion rate M*. Therefore, the inner initial gas disk profile is 
truncated in order to obtain an accretion rate lower that a con- 
stant value of order 3 x 10~ 7 M Q /yr. This allows us to speed up 
the calculation of the disk evolution. 

The initial solid surface density is given by 2d = 
/ D /G/ R/i£o(a/q r 3/2 (Hayashi I198U Weidenschilling et al. 
1997 ) where /d/g is the dust-to-gas ratio of the disk, and /r/i 
is a factor describing the degree of condensation of ices. Its 
value is set to 1 /4 in the regions of the disk for which the ini- 
tial mid-plane temperature exceeds the sublimation of water ice 
(?mid > 170K), and 1 otherwise. The semimajor axes a; ce where 
this temperature is reached as a function of initial gas surface 
density Sq is plotted in fig. [T] For a minimum mass solar neb- 
ula (MMSN) like S (100-200 g/cm 2 ), the iceline is as expected 
found between 2 and 4 AU (Hayashi 1981). Note that in the ac- 
tive disk model we use, the effect of stellar irradiation on the 
temperature structure of the disk is not included. 

2.2. Migration rate 

The migration of the protoplanet occurs in two main regimes de- 
pending upon its mass. Low mass planets undergo type I migra- 
tion (Ward 1997; Tanaka et al. 12002b which depends linearly on 
the body's mass. The prevalence of extrasolar planets have led 
to suspect that the actual type I migration rate is probably sig- 
nificantly lower than currently estimated (Menou & Goodman 
2003 ; Nelson & Papaloizou 2004). For this reason, we allow for 
a arbitrary reduction of the type I migration rate as calculated in 
Tanaka et al. ((2002 ) by a constant efficiency factor fi. 

The migration type changes from type I to type II when the 
planet becomes massive enough to open a gap in the disk. We 
assume that this happens when the Hills radius of the planet be- 
comes greater than the density scale height H of the disk (Lin & 
Papaloizou 1 1986b . Planetary masses were the migration regime 
changes can be low with such a thermal criterium only, as found 
also by Papaloizou & Terquem ( 1999) who use a similar condi- 
tion. This is especially the case as due to disk evolution, the disk 
scale height H decreases with time, so that the minimal mass 
needed to open a gap decreases. This effect is emphasized by 
the fact that our disk model does currently not include irradia- 
tion, so that especially towards the end of disk evolution, H gets 
smaller than in a disk including it, and smaller planets can open 
a gap (Edgar et al. 2007 ). The order of magnitude we obtain is 
however consistent with the one derived from Armitage & Rice 
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Fig. 1. Position of the iceline function of the initial gas 

surface density So at 5.2 AU (upper three lines). It corresponds to 
initial r m ;d of 170 K. The iceline is plotted for three values of a: 
0.01 (dashed line), 0.007 (solid line) and 0.001 (dotted line). The 
lower three lines correspond to initial T m id of 1600 K, roughly 
the evaporation temperature of rock. The rockline a roc k is how- 
ever not taken into account in the nominal model, due to the 
difficulty defining its relevant location, as disk evolution is very 
rapid close-in and irradiation effects might be important (cf. pa- 
per II). 



(2005 ), since they give a gap opening condition (including the 
effect of viscosity) of M p \. dnet /M* > a l ^ 2 (H/a ? \ anet ) 2 . In our sim- 
ulation, the transition typically occurs when the aspect ratio of 
the disk has become tiny, between 2 and 3%, meaning a tran- 
sition at tens of Earth masses. We note that Crida et al. (2006) 
have derived a new criterion for gap opening which depends on 
both the disk aspect ratio and the Reynolds number. Using such 
a modified transition mass has some influence on the planetary 
formation tracks (see sect. |5.3.5[ ). 

Type II migration (Ward [T997l ) itself comes in two forms: As 
long as the local disk mass is large compared to the planet's mass 
A^pianet (called "disk dominated" migration in Armitage 2007 ), 
the planet is coupled to the viscous evolution of the disk and its 
migration rate is independent of its mass. The planetary migra- 
tion timescale is then the same as the gas viscous timescale (e.g. 
Ida & Lin 2004a). Once the local disk mass and the planet's mass 
become comparable migration slows down (Lin & Papaloizou 
1986) and eventually stops. Due to the inertia of the planet the 
disk can no longer deliver the amount of angular momentum 
necessary to force the planet to migrate at the gas' radial speed 
(e.g. Trilling et al.[T998 called "planet dominated" migration in 
Armitage l2507l) . 

As Armitage (120071 ) and Thommes et al. (2008), we have 
found that this braking phase plays a key role in determining 
the final semi-major axis of massive planets. The reason for this 
can be seen in fig. [2] where TLa 2 is plotted as a function of time 
and semimajor axis for an example disk evolving under the in- 
fluence of viscosity and photo-evaporation. The quantity 2£a 2 
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Fig. 2. Example of the evolution of 2Sa 2 . Seven different mo- 
ments in time are plotted (time in units of Myrs). At 4.6 Myr the 
disk has nearly vanished, and the calculations are stopped. 



serves as the measure of the local disk mass to which the planet's 
mass is compared (Lin & Papaloizou 1986 ; Syer & Clarke 1995 1 
Armitage EUOTT ). 

The plot shows that except at the very end, 22a 2 always in- 
creases with a. Initially, in the region where most giant planet 
begin their formation (~5-10 AU), a mass of at least a few Jupiter 
masses is needed to get into the braking phase. However, after 
1-2 Myr, which is the typical timescale to build protoplanets that 
have a sufficient mass to migrate in type II, a mass of the order of 
~ 10M ffl at ~ 1 AU is already sufficient to enter into the slower 
planet dominated type II mode. Combined with eq.[3] fig.[2]also 
shows that once the braking starts, it steadily increases as 22a 2 
decreases with decreasing a and with time, while the planet's 
mass can continue to grow. 

A consequence of the temporal evolution of the gaseous disk 
is the slowing down of type II migration rate thus providing a 
natural mechanism to halt planets at intermediate distances. At 
very small distances from the star (< 0.1 AU), other, special 
stopping mechanisms might additionally be at work (Lin et al. 
[1991)1 

As in Alibert et al. ( 12005 ab . the migration rate during the 
braking phase (M planet > 2S(a p i anet )a 2 lanet ) is calculated as 



da 



planet 



3v ^planet. f )ap lanet 



dt 



"planet 



M, 



planet 



(3) 



where a p i anet is the semimajor axis and v the viscosity. This is 
a modification of eq. 64 in Ida & Lin (2004a) in the sense that 
we use as Edgar (2007) the disk properties like X directly at the 
planet's position and not at the radius of maximum viscous cou- 
pling, assuming that for the rather high a we are using, wave 
dissipation and thus angular momentum exchange will occur es- 
sentially in the proximity of the planet (Lin & Papaloizou 1984). 

It can finally also be noted, that a slowing down of da/dt oc 
2fl 2 lanet /M p i anet explains naturally why larger planets should stop 



further out, provided that Ea 2 increases with a. This behavior is 
indicated by observations (Zucker & Mazeh 2002) 



2.3. Protoplanet structure and evolution 

The structure of the forming planetary envelope is calculated by 
solving the standard equations of planet evolution as in Alibert et 
al. (f2005a), but assuming that the luminosity of the envelope L is 
uniform, and equal to the accretion luminosity of planetesimals: 



GM cnre M, 



core J " core 



Rc 



(4) 



In this equation, M core is the accretion rate of planetesimals, 
M core is the mass of the core of the planet, and R cole its ra- 
dius. The accretion rate is calculated according to Greenzweig 
& Lissauer (|1992| l, using the same prescription for the planetes- 
imal random velocities Vdisp as in Pollack et al. (1996). We as- 
sume a M cole which is independent of migration, for the reasons 
given in Ida & Lin (2008). We therefore do not explicitly com- 
pare the timescales of planetesimal random velocity excitation 
with the migration timescale to see whether the protoplanet acts 
as a "predator" or "shepherd" (Tanaka & Ida 1999). It should in 
any case be noted that if the type I efficiency factors fi are under- 
stood as a consequence of a "random walk" type migration (e.g. 
Nelson & Papaloizou 2004), where the single modifications of 
the orbit occur on a short timescale, then we are in a regime that 
remains yet to be explored in details (Daisaka et al. 2006]). In 
the slow, planet dominate type II regime, where planets are mas- 
sive, accretion of planet esimal s is usually no more important, as 



ejection dominates (sect. 5.1.4 1. When calculating M core , we take 
into account the focusing effect of the planetary envelope (Inaba 
& Ikoma 2003 ). The iterative procedure to get the effective cap- 
ture radius of the planet R capt is the same as in Pollack et al. 
(1 19961 1. In these calculations the effects of ablation are included 
(Benvenuto & Brunini 2008 ). We have found that ignoring the 
focussing effect of the envelope leads to a planetary population 
with similar general properties, but where roughly only half as 
many giant planets can form. 

The assumption that L is uniform and due to planetesimal 
accretion only constitutes a major difference between the models 
in this work and the ones of Alibert et al. (2005a): we do not 
take into account the exact location of the energy deposition of 
infalling planetesimals in the envelope, nor the energy released 
by the contraction of the envelope. Tests have shown that these 
assumptions do not strongly affect the formation, as also shown 
by Rice & Armitage ( |200l1 ). 

Solving the structure equations using the local disk temper- 
ature and pressure as external boundary conditions gives us the 
gas accretion rate of the planet as well as the critical mass for gas 
runaway accretion. Note that Miguel & Brunini (2008 ) have re- 
cently shown that the large uncertainties affecting the constants 
used in parameterized gas accretion laws as in Ida & Lin |2004a) 
lead to large variations of the predicted final planetary mass dis- 
tribution. 

The afore-mentioned method is valid as long as the disk can 
supply enough mass to keep the outer radius equal to the Hill (or 
the accretion) radius, i.e. if the gas accretion rate deduced from 
the envelope structure calculation is below the maximum rate at 
which gas can be delivered by the disk onto the planet. 

This latter quantity can be influenced by the response of the 
disk on the planet's tides. Indeed, hydrodynamical simulations 
(Lubow et al. [19991 D'Angelo et al. [2002b have shown that, 
when the planet opens a gap in the disk, the accretion rate of gas 



C. Mordasini et al.: Extrasolar planet population synthesis I 



5 



is highly reduced. However, it has also been shown (see Kley & 
Dirksen 2006) that when the mass of the planet becomes of the 
order of 3 - 5 (M^ is the mass of Jupiter ^ 318 M e ), the 
disk-planet system can undergo a dynamic instability, leading to 
a substantial increase of the accretion rate of gas. For that rea- 
son, we assume in this paper that the planetary gas accretion rate 
in the disk limited case is simply equal to the accretion rate in 
the disk, namely 



dM, 



planet 



dt 



M d 



isk 



?>nvL. 



(5) 



This setting constitutes another difference to the models in 
Alibert et al. (I2005al >. where we had limited the planet's accre- 
tion rate across a gap according to Veras & Armitage d2003b . 

As mentioned in Alibert et al. (2005a ), since we are primarily 
interested in the mass and semi-major axis evolution of forming 
planets, the planet internal structure is no more calculated once 
the limiting accretion rate Mdisk is reached. Therefore, in this 
second phase, we can no more explicitly compute the capture 
radius R capt of the planet. Rather, it is simply assumed to scale 
with the core radius, i.e. the ratio R CR p t /R C or e is kept constant. As 
a consequence, the amount of solids accreted after the limiting 
accretion rate as been reached is uncertain. This affects however 
mainly large, gas dominated mass planets, and not low mass ones 
(like "Hot Neptunes"). 

2.4. Limitations of the model 

From its conception, our model is well suited to describe the 
formation of giant gaseous planets, but only to a lesser extend the 
formation of very low mass (terrestrial) planets. This is mainly 
due to three assumptions that are made. Additionally, our model 
is a formation, not an evolutionary model (i 2.4.4 1. 



2.4.1. Initial embryo mass 

First, we always assume an initial seed embryo mass of 
A^emb,o=0.6 M @ , similar to Pollack et al. d 1996b or Bodenheimer 
& Pollack (1 19861 1. as at such a mass the opacity in the enve- 
lope becomes sufficient to justify the diffusion approximation 
for the radiative flux (Bodenheimer & Pollack 1986 ). This as- 
sumption implies that our models are only valid for planets with 
final masses exceeding this value by some significant margin. 
More quantitatively, the assumption of M em b,o=0.6 M e is rea- 
sonable when the local isolation mass M; so (cf. §4.4 1 is signifi- 
cantly larger than 0.6 M®. For disks similar to the MMSN, M iso 
is larger than 0.6 M @ only beyond the iceline aj ce (e.g. Lissauer 
& Stewart [T993l . 



2.4.2. Growth after disk dispersal 

Second, we stop the calculations when the gas disk has disap- 
peared (or when the planet has migrated close to the sun, cf. 
below), and ignore all processes taking place later on. While for 
giant gaseous planets this assumption is reasonable (except for 
effects due to the concurrent growth of many planets, see below), 
inside the iceline, for disks similar to the MMSN, growth from 
M; so ~ 0.01 - 0.1 M e to the final masses occurs through giant 
impacts on timescales that exceed by roughly one order of mag- 
nitude typical gas disk lifetimes (Goldreich et al. 2004a). Thus, 
for terrestrial mass planets, growth after the dispersal of the gas 
disk is of large importance. Planetary accretion proceeds at a 
slower pace at larger distances (i 4.5 i, so that our assumption of 



an essentially completed formation of the planets at the time of 
disk dispersion is also not fulfilled for the formation of gas-free 
ice giants at large distances (a > 10 - 30 AU, Ida & Lin 2004a). 

We therefore caution that our synthetic planetary popula- 
tions are incomplete for masses less than a few earth masses for 
a < fl; ce and less than a few 10 M e for a > a; ce . We stress how- 
ever, that the fact that the vast majority of seed embryos do not 
become giant planets (cf. below) is not an artifact of the model 
but a consequence of the fact that the most common protoplane- 
tary disks allow only the formation of relatively low mass plan- 
ets. 



2.4.3. One embryo per disk approach 

An important limitation of the model that is linked to the last 
point is the fact that we follow, as in previous similarly detailed 
giant planet formation models (Pollack et al. 1 19961 Alibert et al. 
2005a), the growth of only one embryo per disk, revolving on a 
circular orbit. 

In reality it is clear that more than one embryo will emerge in 
the same protoplanetary disk, which typically start out forming 
in rapid succession and close proximity (Thommes et al. 2008). 
This multiplicity can have several kinds of effects during for- 
mation, as shown by Thommes et al. (2008): Gravitational in- 
teractions between forming planets can modify their migration 
rate during formation, in particular by locking into resonances. 
Moreover, similar interactions can lead to modification of the 
semi-major and eccentricity distributions, also after the forma- 
tion process itself (see also Adams & Laughlin 12003 1 . In ad- 
dition, planets forming in the same protoplanetary disk act as 
competitors for planetesimals and gas accretion, where e.g. one 
planet can cut off the gas supply of the other ones. The compe- 
tition for planetesimals accretion was in particular addressed by 
Alibert et al. d2005bl > in the case of the Solar System formation, 
and of the HD 69830 system (Alibert et al. 120061 ). In these two 
models, the internal structure of forming planets were calculated 
(as in the present models), but not the gravitational interactions 
between forming planets. In a different, and nicely complemen- 
tary approach, Thommes et al. (2008) have calculated a large set 
of multi-planet formation models (following a much larger num- 
ber of embryos compared to the two afore-mentioned studies but 
without determining the internal structure of planets) accounting 
for both competition for gas and solids accretion, and gravita- 
tional interactions between planets. Including several embryos 
while keeping the detailed physics describing one single embryo 
is a difficult, but important step to be taken in future models. 



2.4.4. Planets very close to the star 

Other complications arise when planets migrate very close to 
their host star (< 0.1 AU). This close to the star, the disk struc- 
ture is more complex than further out due for example to mag- 
netic field effects (Lin et al. 1996 ) or tidal interactions (Trilling et 
al. 120021) . which influence the formation of planets entering this 
zone, by altering the accretion or the migration rate (Papaloizou 
& Terquem 11999b . Also after formation, very close-in planets 
can be subject to mass loss by evaporation (Vidal-Madjar et al. 
120031 Baraffe j2004). Planets in great proximity to the star can 
thus only be described if a more detailed disk model is used, if 
additional stopping mechanisms for migration are taken into ac- 
count and finally, if a subsequent evolutionary model for evapo- 
ration (Baraffe et al. 2006) is included. These complications also 
illustrate the importance of discovering more (small) planets at 
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"safe" distances from the star, as they are a more direct constraint 
on formation models. 

Our model currently doesn't include any of the afore- 
mentioned effects. We therefore simply stop the calculations 
when a planet of mass M p i anet has migrated to a t0 uch, which 
is defined as the semimajor axis where the inner boundary of 
the planet's feeding zone touches the inner boundary of our 
computational disk at a m j n =0. 1 AU, ("the feeding limit") i.e. at 
atouch = fl m in/(l - 4(M p i anet /(3M t )) 1/3 ). If a planet has migrated 
to fltouch. all we can state is that its final semimajor axis would be 
< a t0U ch (it is also possible that it eventually would have fallen 
into the host star), and what its mass at a to uch was. 

3. Monte Carlo method 

The basic idea of using a Monte Carlo method to synthesize 
planetary populations is to sample all possible combinations 
of initial conditions (protoplanetary disk mass, metallicity, etc.) 
with a realistic probability of occurrence. This leads to all pos- 
sible final outcomes of the formation process (i.e. planets) also 
occurring with their relative probabilities. We first explain the 
general six steps procedure that we used. 

In the first step, we have identified four crucial initial con- 
ditions, and studied the domain of possible values they can take 
Some other initial conditions had to be kept constant dur- 



(?3.1 



ing the synthesis of one population, for simplicity or computa- 
tional time restrictions ({j |3.2| i. In the second step, we have de- 
rived probability distributions for each of the four Monte Carlo 
variables (Q. In the third step, we draw in a Monte Carlo fash- 
ion large numbers of sets of initial conditions. The forth step 
consists in using the formation model for each set of initial con- 
ditions, giving the temporal evolution of the planet (formation 
tracks, §5.1[ ) as well as its final properties (mass, semimajor axis, 
composition etc., 5 5.2 1. 



Many of these synthetic planets would remain undetected by 
current observational techniques. So, to be able to compare the 
synthetic planet population with the observed one, we apply in 
the fifth step a detailed synthetic detection bias (paper II). In this 
way, we obtain the sub-population of observable synthetic plan- 
ets. Ultimately, in the sixth step, we have performed quantitative 
statistical tests (paper II) to compare the properties of this ob- 
servable synthetic exoplanet sub-population with a comparison 
sample of real extrasolar planets. 

3.1. Monte Carlo Variables 

We use four Monte Carlo variables to describe the varying initial 
conditions for the planetary formation process. Three describe 
the protoplanetary disk and one the seed embryo. 

1 . The dust-to-gas ratio in the protoplanetary disk f D /c deter- 
mines (together with So) the solid surface density. Models with 
/d/g between 0.013 and 0.13 were computed. Combined with 
the domain of Eo> this corresponds to initial solid surface den- 
sities at ao = 5.2 AU between 0.65 and 130 g/cm 2 . For com- 
parison, the MMSN has a value of approximately 2.5 g/cm 2 
(HayasMQMD. 

2. The initial gas surface density So at 5.2 AU gives the 
amount of gas available. Values between between 50 and 1000 
g/cm 2 were used. The MMSN is estimated to have had a value 
of about 100-200 g/cm 2 (Havashi fTWTI 

3. The last variable that characterizes a disk is the rate 
at which it looses mass due to photoevaporation M w . For the 
population presented below, it was allowed to vary between 
5 X 10- 10 M o /yr and 3 x 10" 8 M o /yr. 



4. Finally, the initial semimajor axis of the seed embryo 
within the disk, a s t art , is the fourth variable. It can take values 

Of 0.1 < fl s tart <20 AU. 

3.2. Parameters 

Some other initial conditions of the model were kept constant 
for all planets of given population. We mention only the most 
important parameters here. More details can be found in Alibert 
et al. d2005ab . For the nominal population discussed in ^5] we 
use a viscosity parameter a for the disk model of 0.007 and an 
efficiency factor for type I migration fj of 0.001. The influence 



of these two important parameters is briefly discussed in 5 5.3.3 



and will be further considered in forthcoming publications. In 
this and the companion paper the mass of the central star M, is 
kept constant at 1 M Q . 

4. Probability distributions 

In the next step we determine the probability of occurrence of 
a certain combination of initial conditions. In the ideal case, the 
probability distributions for all our variables would be derived 
directly from observations. Unfortunately, in reality, this is not 
possible either because in some cases observations do not exist 
or, even if they exist, a certain amount of modeling is necessary 
to to extract the distributions from the observations. 



4.1. Dust to gas ratio / D /g - [Fe/H] 

To establish a link between the dust-to-gas ratio /d/g which is 
the computational variable required by our model and the corre- 
sponding observable, the stellar metallicity [Fe/H], we assume: 
(1) the stellar content in heavy element is a good measure of the 
overall abundance of heavy elements in the disk during forma- 
tion time. Support for this assumption comes from the small dif- 
ferences between solar photospheric and meteoritic abundances 
(Asplund et al. 2005), (2) a scaled solar composition and (3) a 
negligibly small influence of the change of the relative heavy 
element content on the relative hydrogen content in the compar- 
atively small [Fe/H] domain of interest for planet formation in 
the solar neighborhood (-0.5 <[Fe/H]< 0.5). Then, similar to 
Murray et al. (1200 11 1, we can write 



/d/g 
/d/g,g 



(6) 



where /d/g,o ls the dust to gas ratio corresponding to [Fe/H]=0. 
This formula implies that we assume that iron is a good tracer for 
the relevant overall amount of solids available for planet forma- 
tion. Robinson et al. (2006) have found that at given iron abun- 
dance, planet host stars are enriched in silicon and nickel over 
stars without planets, indicating that the above relation is a sim- 
plification. 

Measurements of the heavy element abundance in the Sun 
yield the amount (for complete condensation) of high Z material 
that existed initially in the form of uniformly mixed fine dust 
grains. However, what is relevant for our simulations is the con- 
centration of solids in the innermost 20 AU of the disk at a later 
stage, namely when the dust has evolved into the 100 km plan- 
etesimals used in our model. 

As has been shown by Kornet et al. (2001), the transition 
from the very early dust phase to the later planetesimal phase 
involves a number of coupled mechanism of dust-dust and dust- 
gas interactions like dust settling to the midplane, dust growth by 



C. Mordasini et al.: Extrasolar planet population synthesis I 



7 



Table 1. Parameters p and cr describing the Gaussian distribu- 
tion of [Fe/H], for different observation samples. FV05 stands 
for Fischer & Valenti d2"005l 



Source 




a 


Nordstrom et al. 120041 


-0.14 


0.19 


Fit to CORALIE planet search sample 


-0.02 


0.22 


Fit to FV05 planet search sample 


0.05 


0.21 


Fit to FV05 volume limited sub-sample 


-0.05 


0.26 



coagulation and radial drift. This leads to a redistribution of the 
solids within the disk, which can in turn have important effects 
on planetary formation (Kornet et al. 2005). The key point is 
(Kornet et al. 2004) that these processes lead to an increase of 
the solid to gas ratio in the inner (< 10-20 AU) planet forming 
regions of the disk by advection of solids from the outer parts 
where a lot of mass resides. The factor of increase from the initial 
"(dust-)"/D/G to the "(planetesimal-)"/D/G varies depending on 
the initial conditions and is not completely uniform across the 
inner disk, but Kornet et al. (2004) typically find values of 2 to 
4. 

To take this effect into account at least to first order, we set 
the effective planetesimal f D /c in the inner disk to a value about 
3 times higher than the value that is inferred from the solar pho- 
tosphere. Unfortunately, the latter one has been debated recently: 
Anders & Grevesse ( 1989 ) found a value for the protosolar Z of 
0.0189, whereas a more recent study (Lodders 2003 ) indicate 
lower values (Z=0.0149). Multiplying this value by roughly 3, 
gives a Jd/g,o of ~ 0.04, the value that we regard as nominal in 
our simulations. Note that if dust redistribution mechanisms are 
at work with similar consequences in all disks, we can still use 
the observed [Fe/H] distribution to scale / D /g for other stars. 

Next, we determine the probability distribution for [Fe/H]. 
The distribution of the metallicity of solar like stars in the solar 
neighborhood has been well studied (e.g. Nordstrom et al. 2004), 
and can be well approximated by a Gaussian distribution (mean 
fi and dispersion cr), i.e. 



p([Fe/H]) = 



1 



(7) 



The probability density for /d/g is then found by using eq. 
[6] and the fundamental transformation law of probabilities, 
\ P ([Fe/H])d[Fe/H]\ = \p(f D/G )df D/G \ leading to 



pifo/c) 



log(e) 



(i°g<-fo/G Hjr 



(8) 



JD/GO" V27T 

: log(/b/G,o) + P- This is simply a lognormal distribu- 



where ft 
tion. 

Since we aim to quantitatively compare observations and 
theoretical calculations and compare, for example, real and the- 
oretical detection probabilities, we need the metallicity distri- 
bution of a sample of stars that are actually included in planet 
searches. Such samples can have a different metallicity distri- 
bution than volume limited samples as shown by Fischer & 
Valenti (2005 ) who find that their planet search sample is shifted 
by ~ 0.1 dex towards higher metallicities relative to a volume 
limited subset. The CORALIE planet search sample (Udry et 
al. 12000) consists of about 1650 G and K dwarf within 50 pc, 
and is volume limited. In fig. [3] the [Fe/H] distribution (ob- 
tained from calibrations to the CORALIE cross-correlation func- 
tion) for about 1000 non-binary, slow rotating stars within the 
CORALIE search sample is plotted (Santos et. 2003 ). 



in 
fi 
<D 
Q 

o 



1.5 



1 - 



0.5 




■0.5 
[Fe/H] or [Me/H] 



0.5 



Fig. 3. Solid lines: Histogram and Gaussian fit to the [Fe/H] dis- 
tribution of the CORALIE planet search sample. Dotted line: For 
comparison, distribution of photometric metallicities [Me/H] of 
Nordstrom et al. {2004). Probability densities are given, i.e. the 
area under the curves has been normalized to unity in each case. 



To derive an analytical [Fe/H] probability distribution, we 
have assumed that the CORALIE data is Gaussian, and per- 
formed a non-linear least square fit to obtain the parameters de- 
scribing this distribution, p and cr. The results are p = -0.02 « 
0.0 and cr = 0.22 (Table [TJ which we have used in the cal- 
culations. The fit is also plotted in fig. [3] It approximates the 
data very well between -0.45< [Fe/H]<0.4. We have repeated the 
same procedure also for the planet search sample, and the vol- 
ume limited sub-sample described in Fischer & Valenti (2005 ). 
The results are also given in table [T For these two data sets, the 
Gaussian fit is however somewhat ess good, especially for the 
volume limited sample. 

The distribution derived from the CORALIE survey com- 
bined with our fiducial value for /d/g,o> ar, d the range of 
7d/g> result in the following range of stellar metallicities: 
-0.49 <[Fe/H]< 0.51. This range includes almost all known 
planet hosting main sequence stars in the solar neighborhood. 
Finally, we note that despite the complications of linking /d/g to 
[Fe/H] mentioned here, one can reasonably assume that it is the 
variable that is most directly constrained by observational data. 

4.2. Gas surface density S - Disk Mass M disk 

The probability distribution of gas surface densities So can also 
be at least partially inferred from observations of star forming 
regions. The flux density of thermal continuum emission orig- 
inating from cold dust orbiting young stellar objects (YSO) at 
millimeter and submillimeter wavelengths allows an estimate of 
the total dust disk mass (Beckwith et al. 1 1990b . By adding a gas 
content typical for the interstellar medium (Beckwith & Sargent 
1996; Andrews & Williams 2005), total disk masses M^sk are 
found. The study of many YSO forming concurrently in a given 
star formation region then gives distributions for Mdisk- 
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Fig. 4. Histograms and fits to circumstellar disk masses Mdisk for 
Ophiuchus (solid lines) and Taurus-Auriga (dotted lines). The 
actual gas masses in the computational domain out to 30 AU are 
between 0.004 and 0.09 M . 

Table 2. Parameters describing the lognormal distribution of cir- 
cumstellar disk masses fi and cr, from different sources. The 
value given by Robinson et al. (2006) is obtained by extrapo- 
lating the observed yu = -2.31 of Andrews & Williams (2005) 
back to the assumed initial state, and cr is reduced to 0.25 from 
observed value 0.5. The values of Ida & Lin (2004a) are calcu- 
lated using the value for the MMSN from Hayayshi dl981| l. 



Source 




A*disk(/i)[A4] 


a 


Fit to Taurus 


-1.66 


0.022 


0.74 


Fit to Ophiuchus 


-1.38 


0.042 


0.49 


Robinson et al. (12006b 


-1.3 


0.05 


0.25 


Ida & Lin d2004ab 


-1.48 


0.033 


1.0 



While the details of these distributions differ across known 
star formation regions, most have roughly the shape of a lognor- 
mal distribution for M^, i.e. log(Mdi s k) is Gaussian distributed 
with a mean fi and a standard deviation cr (Andrews & Williams 
2005 ; Robinson et al. 2006). To compute the parameters ju and 
cr we use Beckwith & Sargent (1996 ), who give in their figure 
5 histograms for the distribution of total disk gas masses in the 
Taurus-Auriga and Ophiuchus star formation region, and per- 
form as for [Fe/H] a nonlinear least square fit to obtain the two 
sets of parameters. 

The observational data, as well as the fits are plotted in fig. 
|4j while the values for // and cr are given in table [2 For both star 
formation regions, the disk mass distribution can be reasonably 
well fitted by Gaussians. The resulting distributions are however 
very broad which leads to significant probabilities for disks with 
masses in excess of 0. 1M Q (or even 0.3 M ). For typical thermo- 
dynamical conditions, such disks are self-gravitationally unsta- 
ble in the outer regions (e.g. Ida & Lin 2004a). It is likely that 
this result is due to the presence of the remains of the envelope 



from which the star formed, contributing to the observed flux 
(Andrews & Williams l2005l 

The large spread in masses is also due to the spread in age 
of the observed YSOs in one cluster (Robinson et al. 2006 ) and 
thus an evolutionary effect. Possibly, the spread is additionally 
enhanced by the spread of stellar masses, if the the mass of stars 
and disks are correlated. 

Since our planet formation model includes disk evolution, 
the mass distribution we require is one at a moment in time as 
early as possible, before strong evolutionary effects have taken 
place. Therefore the distribution of Ophiuchus, which is about 
2-3 times younger than Taurus-Auriga (White & Hillenbrand 
2004), is more appropriate for our goal and used as nominal dis- 
tribution for the simulations. 

When converting the observed quantity (disk masses), to the 
corresponding numerical quantity in the model, (the initial gas 
surface density So), we are confronted with a situation similar 
to the [Fe/H] to /b/G conversion (Matsuo et al. 2007): Observed 
protoplanetary disks have physical rad ii a p hy s of typically a few 
hundred AU (e.g. Beckwith & Sargent 1996). For computational 
reasons we however simulate only the inner part of the disk out 
to flmax = 30 AU. Putting the observed disk masses into such a 
small disk would lead to unrealistically high surface densities. 
Therefore, when converting the observed disk masses to So, we 
assume that the mass is contained in a disk with a physical radius 
aphys = 300 AU of which we simulate only the innermost 30 
AU. The outcomes of our simulations are not very sensitive to 
the exact value of the assumed physical outer radius, as the disk 
mass scales only with ya p h ys . Thus, the probability density fol- 
io is 



K£o) 



log(e) 
i=t 

SqCT V27T 



[luglSpHi)- 



(9) 



where Ji-yi- log (4na^ 2 ( ya p h ys - V a min )), which is of a log- 
normal type (Bronstein et al. 1999 ). 

With the boundaries of the computational disk at a m ; n = 0.1 
AU and a max = 30 AU and given the range of So of 50-1000 
g/cm 2 , we cover a range of disk gas masses in the computational 
domain between 0.004 and 0.09 M Q . These masses should be 
stable against self-gravitional collapse (Mayer et al. 2004). 

4.3. Photo-evaporation rate M w - Disk lifetime t diiik 

The values for the photo-evaporation rate M w determines, to- 
gether with the viscosity parameter a, the timescale fdisk on 
which the gas disk disappears. Constraints on its possible val- 
ues can therefore also be inferred from observation. 

Haisch et al. (2001) have shown by near-infrared observa- 
tions of hot dust that the circumstellar disk fraction in young 
stellar clusters is a roughly linearly decreasing function of age 
with essentially all disk having disappeared after ~6 Myr. They 
argue that this timescale should be a tracer of the evolution of the 
bulk of the disk material as well, especially also of the gas disk. 
A very similar time dependence is also found by Hillenbrand 
(120051 

Assuming a uniform distribution in log for M w , we have de- 
termined the bounds of the M w distribution by requiring that the 
lifetime of our synthetic disks for a fixed value of the viscosity 
parameter a and for the Ophiuchus distribution of initial disk 
masses, follows the observed disk lifetime distribution found by 
Haisch et al. (2001 ). To determine the bounds, we have first gen- 
erated a distribution of disk masses as described in sect. 14.21 
We have then set some trial boundaries for the distribution of 
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Fig. 5. Fraction of stars in the model possessing a gaseous disk 
as a function of time for a uniform distribution in log of M w 
between 5 x 1(T 10 and 3 x 1(T 8 M Q /yr (thin solid line). The 
thick solid line is the fit of Haisch et al. (|2001 ) to the observed 
JHKL excess / disk fraction as a function of mean cluster age. 
The error bar at the top right is the overall systematic uncertainty 
in age of these observations, also from Haisch et al. (120011 1. 
For comparison, the fraction is also given if two different dis- 
tributions of M w are used in the model: Uniformly in log be- 
tween 1 x - 1 x 10~ 7 M G /yr (dash-dotted line), leading to 
disk lifetimes that are clearly too short, or uniform in log in 
1 x 10~ 10 - 1 x 10~ 9 M o /yr, leading to disk lifetimes that are 
too long (dotted line). 



M w and drawn values uniformly in log inside these limits. We 
can then calculate the resulting lifetime fdisk^, M v ) for each 
synthetic disk, so that we get a distribution of synthetic fdisk- The 
bounds were then adjusted in an iterative fashion until the ob- 
served and the synthetic disk lifetime distribution are similar. 

It is found that a reasonable fit is obtained by allowing a 
range between 5 x 10" 10 and 3 x 10~ 8 M Q /yr for a = 7 x 10" 3 . 
For other a, the same iterative procedure was repeated. As one 
expects, the higher a, the lower the necessary M w boundaries in 
order to reproduce Haisch et al. (2001). We have for example 
found that for a = 10~ 2 , bounds equal 1 x 10" 10 and 1.5 x 10~ 8 
M G /yr are appropriate. Such values are compatible to the ones 
found elsewhere (Armitage et al. 12003b . 

In fig. [5] the fraction of stars in the model with remaining 
disks using the nominal M w distribution and two comparison 
cases are plotted, together with the fit of Haisch et al. ( 1200 11 1. 
The model and the observed distribution have a similar decrease 
of the remaining disk fraction if the systematic errors of the ages 
in the observational data are taken into account. 

It is clear that the distribution of M w could in reality be much 
more complicated and that also the assumption of a temporally 
constant UV flux might not be justified. One should for example, 
when external UV sources are considered, take into account the 
stellar environment as O or B stars in rich clusters can greatly 



enhance the far-ultraviolet flux (Adams et al. [2004 ; Armitage 
120001) . 

The boundary conditions of our disk model ( §2.1| i are simi- 
lar to the ones presented in Alibert et al. ( 12005 al) . In particular, 
at the outer disk boundary, we assume no mass influx from the 
outer parts of the disk. This assumption could lead to an under- 
estimation of the disk dispersal time, since the outer parts of the 
disk could act as a mass reservoir. However, as explained above, 
we adjust the photoevaporation rate in order to obtain disk life- 
time similar, by construction, to observed ones. With the initial 
power law surface density, the disk accretion rate is found to de- 
creased towards the outer parts of the disk, reaching values be- 
low 10~ 9 M G /yr, which is of the order of, or lower than, the pho- 
toevaporation rate. Therefore, a large part of the mass present in 
the outer disk parts is photoevaporated rather than accreted into 
the inner parts. This means that if we would include inflow, we 
would have to increase the required photoevaporation rates to 
obtain disk lifetimes consistent with observation, but otherwise 
our results would remain similar. 



4.4. Embryo start position a start 

The starting position of the planetary embryos a start inside the 
disk cannot be constrained by observations. Therefore, we derive 
a probability distribution using theoretical arguments only. 

The derivation can be made relatively easily if one assumes 
that during the early phases of planetary accretion radial mo- 
tions can be neglected. In this case, it follows (e.g. Lissauer & 
Stewart fl993t from the restricted 3-body problem that an em- 
bryo can accrete background planetesimals only within its feed- 
ing zone which has a half width of Bi times (Bl » 3 - 5) 
its Hill sphere radius R^. The Hill sphere radius of the plane- 
tary embryo at a semimajor axis a and a mass M em b,o is given 
by R H = (M em b,o/(3M t O)) I/3 a, which is, for fixed M emb ,o, <* 
a. Thus, as confirmed by many different numerical simulations 
(e.g. Weidenschilling et al. 1 19971 Kokubo & Ida 2000) runaway 
bodies emerge with relative separations of their semimajor axis 
A proportional to their semimajor axes a, so A/a = const. This 
situation prevails also later during oligarchic growth stages, al- 
though the numerical value of Bi increases (Ida & Lin 2004a). 

The probability p(a) of starting in an interval da is inversely 
proportional to the spacing A, so 



da da .... 
p(a)da oc — oc — = fllog(a) oc const., 
A a 



(10) 



which means that the probability distribution for the starting po- 
sition is uniform in log, as assumed also by Ida & Lin (2004a). 

While the distribution of the initial position is clear, there re- 
mains the issue that we start our simulations with a seed embryo 
of mass M em b,o - 0.6 M @ . In order to be self-consistent, we have 
to require that this amount of mass is available at the starting 
position. In other word, we have to ensure that for a set of initial 
values (/d/g> an d fl start) there is indeed enough mass avail- 
able to build-up such a seed. For this, we consider the amount of 
heavy elements in the planet's feeding zone 27ra st . drt 2BLRii"Lu to 
calculate the isolation mass M; so (Lissauer & Stewart [19931 ): 



M iso 



£d) 3/2 



(3M„) 1/2 

For the disk profile we use this becomes 

_ (4nB L f D/c f R/ ^ al /2 a l J^/ 2 
(3M«)i/2 



(11) 



(12) 
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Fig. 6. Probability distribution for a s tart- There is a marked peak 
around 3 AU. Outside this distance, the probability distribution 
falls off as expected for a distribution that is uniform in log. The 
fraction inside about 2 AU is due to disks with the ability to form 
seeds sufficiently massive (> M em b.o) inside the iceline as well. 



where / R /i is itself a function of Eq (and a, fi g.[T] i. 

If the isolation mass calculated from eq. [l2|with Bi — 2 V3 
(Lissauer & Stewart 1993J for a set of /b/G, So and a start exceeds 
0.6 M & , we conclude that these are self-consistent initial condi- 
tions for the formation model. 

As M; so oc fl^, the isolation mass criterion sets an inner 
boarder for the possible domain of a s tart- We also require that the 
starting time of the embryo r start is smaller than the disk lifetime 
^disk (cf. the next section jj |4.5| ). In practice, we first calculate for a 
given / D/G , So and M w the range of a start in which M iso > M emb , 
and fdisk > fstait, and then draw a start from that range. 

In fig. [6] the resulting distribution of a start is plotted. Its par- 
ticular shape is the results of the combination of the probabil- 
ity distribution for a start (uniform in log) and the criteria that 
fdisk ^ fstart an d Mso - A^emb.o- The fact that the latter is fulfilled 
only beyond the iceline for most disks, and the property of the 
uniform-in-log distribution to emphasize the smallest possible 
values lead to the peak near 3 AU. Further out, the probability 
decreases, a consequence of the uniform-in-log distribution and 
the increasing f start . Only disks with a high surface density (i.e. 
a concurrently high So an d /d/g) which can form sufficiently 
massive seeds (> M em b,o) inside the iceline contribute to the dis- 
tribution inside ~ 2 AU. The transition to higher probabilities is 
not sharp as the iceline is itself a function of So (fig.[TJi. Note that 
via the M\ so and f^k criteria, the shape of the a start distribution 
depends on the /d/g, So and M w distributions. 

4.5. Embryo start time f start 

Depending upon initial conditions, the time needed to form just 
the initial seed embryo of 0.6 M @ is not short compared to the 
overall disk lifetime. We take this effect into account by intro- 
ducing a time delay f start between the beginning of the disk evolu- 



tion and the time at which put the embryo into the disk. This time 
delay is calculated in a deterministic way as a function of /d/g, 
So and a s tart- It is therefore not an independent Monte Carlo vari- 
able. Similarly to a s tart, the time delay cannot be constrained by 
observation. We therefore again rely on theoretical arguments. 

During the early stages of planetary accretion, the process of 
embryo growth proceeds mainly by accumulation of background 
planetesimals and can be well described in the two-body approx- 
imation (Kokubo & Ida 2002 ). The accretion rate of an embryo 
growing from background planetesimals is thus (e.g. Lissauer & 
Stewart [T9931 



dM e 



dt 



1 + 



esc 
<Hsp J 



(13) 



Here, R em h is the radius of the embryo, v esc the escape velocity 
from it and v^sp the typical random velocity of the planetesimals. 

As the embryo grows, the solid surface density of planetesi- 
mals Sd must decrease correspondingly (Thommes et al. 2003): 



</s D 

dt 



(3M.) 



1/3 



dM„ 



67Taj tm ,B L M emh 



!/3 dt 



(14) 



Looking at the dependence of the embryo growth time on semi- 
major axis, it is found (Steward & Ida 2000; Weidenschilling 
& Davis 2001 ) that f start increases significantly with distance to 
the star. The reason is that both the surface density of planetes- 
imals and the frequency of collisions decrease with distance, as 
the pace on which the later occur scales with Q (Kokubo & Ida 
2002). This increase of the growth timescale with distance is 
well reproduced in numerical simulations where a "wave" of em- 
bryo growth is seen to propagate from the inner parts to the outer 
parts of the disk (Weidenschilling & Davis l200Tl Thommes et al. 
120031 . 

Thus, building up a body of fixed mass, in our case the initial 
seed, takes longer for larger a s tart, a trend that holds everywhere 
but at the iceline where S D increases suddenly. To derive the 
corresponding f start , we integrate numerically equations 13 and 
14 until the mass has reached M em b,o = 0.6 M®. As in the for- 
mation model, we use the procedure of Pollack et al. ( 119961 ) to 
calculate the planetesimal random velocity v<ji S p and assume a 
constant planetesimal size of 100 km. As initial values at t = 0, 
we use a mass of 0.66ra 3 , /5 M 2/5 (Chambers 2006). 

pla iso v 1 ' 

Compared to the formation model which calculates the evo- 
lution of the protoplanet starting with M em b,o - 0.6 M m , we ne- 
glect several things in the calculation of f start : First, migration 
is neglected. Tests have shown that with the type I migration 
efficiency factor / we use for the nominal case, migration of 
bodies this small is indeed negligible. Second, we use just the 
simple particle in a box law for the gravitational focussing from 
eq. 13 instead of the full three body results form Greenzweig & 
Lissauer ( 1992 ), and third, we neglect the increase of the capture 
radius R capt due to the gaseous envelope. The latter assumption is 
appropriate given the small seed masses and the assumed plan- 
etesimal size as shown by our own tests but also by Kornet & 
Wolf (2006 ) or Fortier et al. ( EOOTi 

Fig. [7] illustrates the result of the numerical integrations to 
obtain f start . For a solid surface density at 1 AU of 7 g/cm 2 (~ 
MMSN) and for 35 g/cm 2 (~ 5 x MMSN) four snapshots in time 
of the embryo mass as a function of semimajor axis are plotted. 

To allow comparisons with other models (Thommes et al. 
120031 Ida & Lin l2004al Chambers 120061 we have continued the 
integration up to M; so for this plot instead of stopping at M em b,o 
as we do when generating the initial conditions. The dotted line 
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Fig. 7. Snapshots of the embryo mass (solid line) as a function of semimajor axis at four moments in time for two different solid 
surface densities. The dashed line is the isolation mass. The dotted line is M em b,o = 0.6 M®. The initial solid surface density at 1 AU 
is 7 g/cm 2 (left panel) and 35 g/cm 2 (right panel). It should be kept in mind that this kind of calculation is just needed to generate 
the start time f start when the embryo is put into the formation model. The real evolution of the solid core for M > 0.6 M @ is in 
general much more complex than plotted here. In this figure, we have continued the calculations up to the isolation mass just to 
allow comparisons with other models. 



in fig.[7]is M; so and shows that for the MMSN, the earliest pos- 
sible time to start is somewhat more than 1 Myr, just outside the 
iceline. Then, the domain of possible start positions grows only 
slowly to larger semimajor axes. After 10 Myr, when the gas 
disks will have disappeared (cf. fig. [5}, the largest possible start 
position is still only 6-7 AU. In contrast, in the 5 x MMSN case, 
M; so is larger than M em b,o inside the iceline too, and the earliest 
embryo can start at around 0.8 AU slightly before 0.1 Myr. At 
later times, embryos can start from all semimajor axes. 

Compared to the results of Thommes et al. (2003), our cal- 
culations show a faster growth of the cores, especially at large 
distances, which is due to the different way of computing plan- 
etesimal eccentricities and inclinations in their model, as illus- 
trated by Fortier et al. (120071) . Compared to Ida & Lin d2004at . 
the results are quite similar, even if core growth proceeds at 
large orbital distances somewhat faster in our model. Compared 
to Chambers (2006) one finds that core growth in our model is 
faster than in his simple equilibrium model, but slower than in 
his complete model that is considerably more complex, includ- 
ing e.g. planetesimal fragmentation. 



As mentioned in 5 4.4 we only start embryos in that part of 



the disk where M; so > M em b,o and fdisk ^ fstart- The latter con- 
dition gives an outer bound for possible starting positions. The 
reasoning behind it is that if one of the numerous planetary seeds 
can form while the disk is still present, it would indeed have 
done so, and that it is a candidate to eventually become a giant 
planet observable today. In other parts of the disk, seed embryos 
also form, but they remain very small during the whole presence 
of the gaseous disk. Thus, we aim at minimizing the negative 
side effects of having only one seed per disk on the population 
of giant planets, but at the same time make our populations in- 
complete at small masses (cf. 4 2.4 1. 



For a significant fraction (~ 28%) of the sets of initial con- 
ditions we draw, one or both of the two aforementioned condi- 
tions cannot be fulfilled anywhere in the disk, namely when f D /c 
and/or Eo come from the low tail of their distributions, while M w 
is high. In such cases, no calculations were made, but we keep 
the record of the corresponding initial conditions where the for- 
mation of sizable planets is not possible and correct for them 
when calculating for example overall detection probabilities (pa- 
per II). 



5. Results 

Once all Monte Carlo variables have been drawn, the next step 
consists in computing the formation of the planet correspond- 
ing to these initial conditions. This process can be illustrated by 
means of formation tracks in the mass-distance plane. Except 
where otherwise stated, all results are obtained for a population 
with a =0.007 and fi - 0.001. The reason for this choice is 
that the resulting sub-population of observable synthetic planets 
reasonably well reproduces the observed population (paper II). 



5.1. Planetary formation tracks 

Figure [8] shows formation tracks of about 1500 randomly cho- 
sen synthetic planets. The tracks lead from the initial position at 
a(t = 0) = a SIart and the fixed M(t — 0) = M em ^o to the final po- 
sition marked by a large black symbol when planet growth and 
migration stops. The color of the track indicates the migration 
mode: Red for type I migration, blue for ordinary (disk domi- 
nated) type II migration and green for the braking phase. In this 
phase, planet dominated type II migration occurs (eq. [3j and the 
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Fig. 8. Planetary formation tracks in the mass-distance plane. The large black symbols show the final position of a planet. The shape 
of the symbols is explained in the text. Planets reaching the feeding limit at a to uch (indicated by the long dashed line) have arbitrarily 
been set to 0.1 AU. The short dashed lines have a slope of -n (discussion in jj |5.1.3 1. Each track is color-coded according to the 
migration mode, and small black dots are plotted on the tracks all 0.2 Myr to indicate the temporal evolution of a planet. 



planetary gas accretion rate is given by the rate at which the disk 
can supply gas (eq.[5]l. 

Even if the tracks show a great diversity, one can distinguish 
groups of planets with similar tracks. These groups are due to 
different formation stages that planets might undergo. In the next 
sections, we study representative tracks of four such groups. 

5.1 .1 . Tracks of "Failed cores" 

During the first stage of formation at low masses type I migration 
(red) occurs. Since for this example population type I migration 
is very slow (/i = 0.001), the tracks are almost vertical. Planets 
that have migrated in type I only are represented by filled circles 
infig.|8] 

For most embryos, this first stage is also the final one. Their 
evolution stops at small masses because most initial conditions 
do not allow the formation of more massive planets during the 
lifetime of the disk. Therefore, most seeds (50-75 %, see paper 
II) contribute to building up a large population of "failed cores" 
with M ~ 1-10 M m which, for a point of view of giant planet 
formation, failed to accrete a significant amount of gas. Also 
the population synthesis calculations of Ida & Lin (2004a, 2008 ) 
contain a large sub-population of low mass planets. This is also 



compatible with the non-detection of giant planets around 90 to 
95% of nearby solar like stars. 

In fig. 9] left panel, exemplary formation tracks for a number 



12 for Mu 



of such planets are plotted. As expected from eq. 
"failed cores" can reach larger masses at larger distances. The 
right panel of fig. [9] shows the temporal evolution of the mass 
and semimajor axis of one typical "failed core". This seed starts 
at a start = 3.7 AU in a disk with / D/G = 0.028 ([Fe/H]=-0.15) and 
So = 165 g/cm 2 . This initial position is situated not far outside 
the iceline. For such a solid surface density, forming the initial 
seed takes a significant amount of time (cf. fig.|7]i, namely about 
1.1 Myr. 

As is shown in the right panel of fig. [9] the core then quickly 
accretes all planetesimals in its reach. Gas accretion is of neg- 
ligible importance. At about 1.2 Myr, the mass of the core ap- 
proaches the local isolation mass^j For the remaining 0.2 Myr 
of evolution, the core grows only very slowly. One notes that 
the envelope now becomes somewhat more massive, due to the 
reduced luminosity of the core. Clearly, the evolution of this 



2 From eq.[TT|one would calculate a M iso of about 2.8M e , using B L = 
4. However, as we do not reduce the initial solid surface density by the 
amount of material already in the initial seed, a value larger for the mass 
by about 3/2 x Af em b,o. is obtained. 
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Fig. 9. Planetary formation tracks in the mass-distance plane for "failed cores" (left panel). The black points represent the final 
position of the planets. The thick line starting at 3.7 AU is the track of one prototypical example, for which the right panels shows 
the temporal evolution. Its final position is represented by a large square. In the right panel, the total mass M (solid line), the mass 
of accreted solids Mz (dashed line) and the mass of the envelope M em multiplied by a factor 10 for better visibility (dotted line) are 
plotted (scale on the left) as a function of time t. The temporal evolution of the planet's semimajor axis a is also plotted (dash-dotted 
line, scale on the right). 



planet simply corresponds to the two first phases described by 
Pollack et al. (1996), with the difference that further evolution 
is inhibited by the dispersion of the protoplanetary nebula after 
1 .45 Myr. At this time, we are left with a "failed core", consist- 
ing of about 3.6 M @ of heavy elements, and ~ 0.1 M e of gas. 
The extend over which migration occurred is tiny because of 
fi = 0.001, roughly 0.004 AU, much less than the extent of the 
planet's Hills radius. The fact that further growth is simply in- 
hibited by the disappearance of the gaseous disk is characteristic 
for this type of planet. 

The vast sub-population of "failed cores" is not identical to 
the final terrestrial planet population, expected to be located in 
a similar a - M region. Rather, they represent an earlier mo- 
ment in evolution. "Failed cores" are formed from one large em- 
bryo accreting small field planetesimals while the gas disk is 
still present. Terrestrial planets on the other hand get their fi- 
nal properties from giant impacts between bodies of a similar 
size (several "failed cores") on much longer timescales, a phase 
missing in our model. We expect that after disk dispersal, all the 
"failed cores" of one disk would start to interact gravitationally, 
leading to scattering, ejections and collisions, until the remain- 
ing planets have settled into stable orbits (e.g. Ford & Chiang 
2007; Thommesetal. HDDS). 



5.1.2. Tracks of "horizontal branch" planets 

In some other cases the core grows so large (and does so suffi- 
ciently quickly) that the planet can open a gap in the gas disk 
long before the latter disappears. At this point, the migration 
mode changes to disk dominated type II (blue lines in fig. |HJ, 
which is the second phase. 



After a short transitional phase, planets starting inside about 
4-6 AU begin then to move in disk dominated type II migration 
along nearly, but not completely horizontal tracks at M ~ 7 - 30 
M©. These nearly horizontal tracks are clearly seen in fig. [8] 
forming a "horizontal branch" of planets. We identify the planets 
having had such phase during their formation a posteriori by the 
condition that while in type II migration, one finds d J°^ < 0. 1 . 
Physically this means that migration occurs on a significantly 
shorter timescale than accretion. Planets having passed through 
the "horizontal branch" have their final position marked by tri- 
angles in fig. [8] 

Figure [10] left panel, shows some exemplary formation 
tracks of planets that stay in the "horizontal branch" until the 
gas disk has disappeared (so that they end up at intermediate dis- 
tances), or until they reach the feeding limit (so that they have a 
final position < 0.1 AU). 



The prototypical example of the right panel of fig. 10 is 
formed in a disk whose mass is similar to the mean values 
adopted in this work (Lq = 270 g/cm 2 ). The dust to gas ratio 
is with f D /G = 0.03 somewhat smaller than the mean value. 
This leads to a formation time of the initial seed of just about 

1 Myr. For this disk, the iceline is located at 4.2 AU, therefore 
the planet accretes icy planetesimals at the beginning of its for- 
mation. During the first 1.6 Myr, the formation is similar to the 
one described in classical core accretion papers, namely a rapid 
core formation (up to about 8 M @ , the isolation mass, in 0.15 
Myr) followed by a phase of low mass growth, similar to phase 

2 of Pollack et al. ( 1996). Just before t = 1.7 Myr however, mi- 
gration switches to type II, due to the concurrent growth of the 
planet and the decrease of the disk scale height with time, so that 
also a planet of a relatively low mass can open a gap in the disk. 
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Fig. 10. As fig. [9] but for planets of the "horizontal branch". In this case, the prototypical track (thick line, and large square at the 
final position) starts at a start = 4.7 AU. It ends as a "Hot Neptune" planet in the feeding limit at 0. 1 AU. Its temporal evolution is 
shown in the right panel. 



In the population presented here, we have reduced type I mi- 
gration by a large factor. Therefore changing from the (strongly 
reduced) type I to (normal) type II results in a net increase of the 
migration rate and the planet moves into regions of the disks that 
have not yet been depleted in planetesimals (Alibert et al. 2004). 
This significantly increases the core growth and and hence its lu- 
minosity. The latter translates in a slight decrease of the envelope 
mass at 1.7 Myr. 

Shortly after switching migration type, the planet crosses the 
iceline, which reduces the solid accretion rate (see the little kink 
on the mass lines just after 1.7 Myr). During the remaining ac- 
cretion inside the iceline, another ~ 17 M e of rocky planetesi- 
mals are collected. Comparing this to the total amount of heavy 
elements initially inside the iceline (~ 20 M @ ) shows that the 
planet is quite efficient in emptying the planetesimal disk. As a 
consequence, the final core mass is approximately the sum of the 
mass collected while passing through the "horizontal branch", 
plus the mass of the icy planetesimals accreted during the initial 
in-situ growth phase. 

At roughly 1 .9 Myr, the local disk mass becomes comparable 
to the planet's mass so that the migration rate starts to slow down 
(fig. [2]) and the planets starts to accrete gas. However, migration 
continues (albeit at a reduced rate) and the planet enters at t — 
2.2 Myr into the feeding limit at roughly 0.1 AU, where we stop 
the calculations. 

Contrary to the two remaining types of representative tracks 
described below, the gas accretion rate of "horizontal branch" 
planets is governed during their entire formation by the planet 
i.e. the maximum accretion rate limited by the disk is never 
reached. Since the gas accretion rate is moderate (due to the 
fact that the core luminosity is most of the time quite large), 
the planet at the end of its formation is made of a large core 
(~ 26 M B ) and a significant, but still much smaller envelope 
(~ 6 M m ). We find that the sub-population of close-in, low mass 
planets (M < 30-40 M e ) is characterized by a ratio M env /M core 



that varies between -0.02 (at M ~ l0M e ) to ~ 0.3 (at M ~ 40 
M ffi ), i.e. these planets have a structure roughly comparable to 
Neptune. This could be different for close-in Neptune mass plan- 
ets formed through giant impacts between initially smaller bod- 
ies (Ida & Lin 2008), especially as the impacts could remove 
their anyway tenuous gaseous envelopes. Transiting Neptune 
mass planets around solar like stars can serve as a diagnostic to 
distinguish these different formation channels. The recently de- 
tected transiting "Hot Neptune" HAT-P-llb (Bakos et al.[2009) 
has a radius compatible with a rock/ice core with a 10% H/He 
envelope, whereas a pure rock/ice planet (as well as a miniature 
gaseous planet) are excluded (Bakos et al. [2009] ). This suggest 
that this planet was formed in a similar way as described here. 

The "horizontal branch" is thus the "conveyor belt" by which 
Neptune-like planets are being transported close to the star. 
These "Hot Neptunes" are a sub-population of planets that high 
precision radial velocity surveys (using e.g. HARPS, see Lovis 
et al. 120061 ) now find in increasing numbers. Note that subse- 
quent evolutionary effects, namely evaporation can sometimes 
significantly modify the structure of these bodies (Baraffe et al. 

ESS. 

5.1 .3. Tracks of "main clump" planets 

The fact that the tracks in the "horizontal branch" are not com- 
pletely horizontal, i.e. that growth in mass continues in this phase 
is important for the further evolution of the third group of plan- 
ets, the planets of the "main clump". 

These are planets with final masses mostly between the mass 
of Saturn and three Jupiter masses, and semimajor axes mainly 
between ~ 0.3 and 2 AU (see fig. 13 1. For these planets, the 



core grows to a size that triggers runaway gas accretion while 
the planet is collecting solids as it passes through the "hori- 
zontal branch". Once runaway is triggered, the gas accretion 
rate increases very rapidly. With its rapidly growing mass, the 
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planet soon exceeds the local disk mass and migration enters the 
planet dominated type II regime. The planet leaves the "horizon- 
tal branch" upwards in mass thereby starting its third phase of 
formation (plotted in green in fig. [8} . 

During this final braking phase, the planets migrate at the 
reduced type II rate (e q. [3) , and accrete gas at a rate given by 
the disk's evolution (eq.|5Jl. This has the interesting consequence 
that if we combine these equations, we find that planets move on 
formation tracks which are in the log(a) - log(m) plane straight 
lines with a slope c/ log M/d log a = -it. This behavior is clearly 
visible in the formation tracks. In fact, we have used the criterion 
| d log M/d log a + n\ < 0.1 to identify the braking phase a poste- 
riori. The evolution along these straight lines slows down in time 
as can be seen by the increasing number of small black ticks on 
the track near the final position of the planet (fig. [8}. This is a 
consequence of the concurrent decrease of the gas accretion rate 
due to disk evolution, and the slowing down of migration. At the 
end, the planet has accreted all the outer disk gas that has not 
been photo-evaporated. 

Examples of planets that undergo such a three staged evolu- 
tion are plotted in figure 1 1 Two specific planets have reached 
the feeding limit near 0.1 AU, illustrating how Pegasi planets 
form, even though our model does not further treat the forma- 
tion of "Hot" planets once they have migrated to a t ouch- 

The prototypical planet of the "main clump" in the right 
panel of this figure is formed in a disk with / D /g=0.02, i.e. 
[Fe/H]=-0.3, and £o = 280 g/cm 2 . The planet starts its forma- 
tion beyond the iceline (a st art = 6, fli ce = 4.3 AU), and the 
beginning of its formation (up to 4 Myr) is similar as for the 
"horizontal branch" planet before: the planets empties its feed- 
ing zone reaching an isolation mass of about 6 M m . Just before 
t — 4 Myr, migration switches to type II. Shortly thereafter at 
4.1 Myr the planet crosses the iceline (see the prominent kink 
on the solid line in fig. [TTJ. Finally, the braking phase starts at 
about 4.4 Myr. Switching from type I to type II migration results 
in an increase of both the solid accretion rate and the core lu- 
minosity with a corresponding loss of envelope. The difference 
to the previous case originates from a longer living disk and the 
fact that the planet started at a larger distance from the central 
star. The solid mass available while migrating through the "hori- 
zontal branch" is therefore larger (it scales with r 2 E D oc r° 5 ) and 
the core reaches a mass large enough (21 M @ ) to trigger a run- 
away accretion of gas around 4.6 Myr, significantly before the 
gas disk disappears (or the planet reaches the feeding limit as 
for the "horizontal branch" case). The gas accretion rate then in- 
creases rapidly, but soon (at 4.7 Myr) reaches the maximum rate 
allowed by the disk which is a rather moderate 1 .6 x 10~ 4 M e /yr. 
The disk limited accretion rate then decreases slowly with time 
as the disk evolves. From 4.7 Myr onwards the planet is in the 
braking phase with the characteristic -it slope in the a - M dia- 
gram until the disk disappears at 6 Myr. At this point, a Saturnian 
planet has formed at 0.9 AU with a total mass of about 130 M m 
and a Mz ~ 25M e . The prototypical planet also illustrated how 
the combination of effects that are per se straightforward (chang- 
ing the migration regime, crossing the iceline, disk limited gas 
accretion rate) can lead to a complex formation history. 

It is interesting to note that the "main clump" region of the 
a-M diagram is somewhat overpopulated in the observed popu- 
lation as well. Examples are (among many others) HD 100777b 
(Naef et al. 120071 or HD 142b (Tinney et al. 120021 . 



5.1 .4. Tracks of "outer group" planets 

For starting positions larger than ~ 4-7 AU, the formation tracks 
are different. In the case a high /d/g and £q and a low M w are 
drawn together, the core growth timescale is small compared to 
the disk depletion timescale even at these large distances. As 
the amount of solid material available is large, embryos can then 
grow to a supercritical mass almost in-situ (nearly vertical tracks 
up to several tens of M @ in fig. [SJ, without the need of collect- 
ing solid material by migration. Significant migration can nev- 
ertheless occur but only in the planet dominated type II mode 
and when the gas accretion rate is regulated by the disk. Such 
planets which do not have a "horizontal branch" phase are rep- 
resented in fig. [8]by filled squares. The final semimajor axes of 
the in-situ supercritical planets are outside ~ 0.4-1 AU, but over- 
lap with "main clump" planets. Sometimes "outer group" planets 
become extremely massive "Super Jupiters", with masses more 
than one order of magnitude larger than that of Jupiter. 



Examples of such tracks are plotted in figure 12 Comparing 
these tracks with the ones of the "main clump" shows that there 
is a continuous transition between the two types. The prototyp- 
ical "outer group" planet in the right panel of fig. 12 is formed 
in a massive, metal-rich (Eo = 500 g/cm 2 , [Fe/H]=0.3) and long 
lived disk. Due to the large a stal -t (10 AU), the isolation mass at 
the initial location is then very large (around 150 M B ). The seed 
starts at 1.1 Myr, and rapidly switches to type II migration at 
1.2 Myr. The critical mass (around 25 M @ ) is attained already 
at 1 .24 Myr at a position very close to a stMt , triggering a rapid 
accretion of gas. At 1.25 Myr, the gas accretion rate reaches the 
value limited by the disk of initially initially 2.5 x 10~ 3 M e /yr. 

Note that the final amount of heavy elements in this planet 
is quite uncertain. Indeed, "outer group" planets undergo a large 
part of their formation history in the disk limited gas accretion 
phase. As mentioned above, the internal structure of planets dur- 
ing this phase is no more calculated, thus their solid accretion 
rate (via /? cap t) is uncertain. As a consequence, the final Mz found 
here (around 250 M @ ) is unsure, and only a lower boundary 
(~ 25 M m ) is actually well determined. Internal structure models 
(Baraffe et al. 2008) do however require comparable amounts of 
heavy elements to reproduced the observed mass-radius relation 
of the transiting "Hot Super Jupiter" HD 147506b (Bakos et al. 
120071) . Thus, even if this planet is in terms of semimajor axis 
very different than the planet discussed here, it still illustrates 
that objects with very large amounts of solids in their interior 
seem to exist. Figure 12 shows that the forming planet is even- 



tually strongly dominated by gas (about 10 Mq^_ are accreted). 
Therefore, the uncertainty on the final total mass (which is of 
primary interest in this study) is anyway low. 

The figure also shows that M z (with all the uncertainties dis- 
cussed above) reaches its final value between 8 and 9 AU al- 
ready. The planet could in principle continue to accrete all the 
planetesimals down to its final location at 4.5 AU. However, the 
planet is so massive at this stage that nearly all planetesimals are 
actually ejected and not accreted (Ida & Lin 2004a). We how- 
ever find that the large size of forming planets makes them less 
effective in ejecting planetesimals than old, compact objects of 
the same mass. 



5.2. Mass-semimajor axis diagram 

The planetary formation tracks illustrate how planets grow in 
mass and migrate to their final position. In this section, we dis- 
cuss some aspects of the distribution of final masses and posi- 
tions. 
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Fig. 11. As fig. 9] but for planets passing through the "horizontal branch" to become members of the "main clump". Here, the 
prototypical embryo (thick line, square at the final position, temporal evolution in the right panel) eventually leads to a Saturnian 
planet situated at 0.9 AU. 
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Fig. 12. As fig. [9) but for in-situ critical cores which end up in the "outer group". The prototypical embryo starts in a massive, metal 
rich, long lived disk which allows the formation of a very massive ~ IIM7, "Super Jupiter" ultimately located at 4.5 AU. 



5.2.1 . Diversity of planets 



Figure 13 shows the mass and the semi-major axis of N t 



synt 



50000 synthetic planets. The first striking result is that core ac- 
cretion allows for a very diverse planet population. The final 
mass of the planets varies between the smallest possible mass 
(0.6 M m ) and a few extremely massive ~ 40 planets. The fi- 
nal position varies between the innermost possible radius (« 0.1 



AU) to about 18 AU. We conclude that the observed diversity 
of extrasolar planets is within the core accretion paradigm a nat- 
ural consequence of the observed diversity of the properties of 
the protoplanetary disks, which we have varied within the ob- 
served range. It is therefore obvious that a better understanding 
of the properties of the protoplanetary disks, especially the inner- 
most region, has important implication for this planet formation 
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Fig. 13. Final mass M versus final distance a of iV synt ~ 50 000 synthetic planets of the nominal planetary population. The feeding 
limit at a t0 uch is plotted as dashed line. Planets migrating into the feeding limit have been put to 0. 1 AU. As a touc h gets very large 
for M > 20M7j_, also a few extremely massive planets are in the feeding limit which should however be regarded as a simulation 
artifact because our simplification of putting planets that reach the feeding limit to 0. 1 AU ceases to be justified. 



theory. In our population synthesis, a number of parameters are 
kept fixed at all times. We can speculate that in nature most of 
these quantities also fluctuate to some degree. The core accre- 
tion mechanism should therefore be able to produce planets of 
an even greater diversity than found here. On the other hand, the 
question of whether the core accretion can explain all the planets 
(cf. Matsuo et al. 120071) is a much more difficult one. The models 
are probably not yet mature enough and observations are still too 
incomplete to allow a definitive statement. 



Figure 13 reveals that the synthetic planets are not randomly 
distributed inside the mass-distance plane. Instead, various con- 
centrations, bars and depleted regions can be distinguished. One 
can also study the form of the quite well defined envelope filled 
by the synthetic population, and why certain regions remain 
empty. When comparing fig. 13 with actual discoveries, one 



should however bear in mind the incompleteness of the model, 
affecting in particular low mass planets ({ 2.4 1. 



5.2.2. The limiting envelope 

At large masses, the planetary population is bounded by the 
largest overall mass for a given semimajor axis, M max (a). 
Outside a ~ 3 AU, M max is a decreasing function of a, falling 
to about 100 M m at ~ 20 AU. The reason is that the time needed 
to build up a massive core at large distances becomes compara- 
ble and ultimately longer than the gas disk dispersion timescale 
(Ida & Lin 2004a ). See sect. |5.3.4| for processes which could 
modify this behavior. 

Inside 3 AU, the behavior of M max is inverted, i.e. M max is an 
increasing function of a. This means that there is an absence of 
very massive planets at small orbital radii. This is what has been 
discovered in the observed extrasolar planet population if only 
planets around single host stars are considered (Udry et al. 2003 ; 
Zucker & Mazeh 2002). With the initial surface density profile 
used here (S oc cT 3 / 2 ), we find that inside 3 AU, M max scales 
approximately as a 3 ^ 4 (as M; so ), provided that type I migration 
is slow, see 5 5.3.3 For populations obtained with higher type I 



migration rates M max is flat inside ~ 3 AU. Future discoveries 
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of a very large number of single giant planets out to several AU 
(Ge et al. 2007 ) around single stars will help to define better the 
exact shape of M max (a). 

5.2.3. Structures in the a - M diagram 

The different phases of planet formation and migration that were 
identified in the formation tracks leave traces also in the final 
a - M of the planets. One can distinguish the "failed cores", the 
"horizontal branch", the "main clump", and the "outer group" 
planets. As a new feature, fig. 13 however also shows a deple- 
tion of planets with masses between 30 to 100 M m . This is the 
analogue of the "planetary desert" first discussed by Ida & Lin 
(2004a). Compared to their results, the depletion is much less 
severe in our simulations. 

The reason for this difference is difficult to pinpoint, as both 
formation models differ in many aspects, but it is at least par- 
tially due to the way the maximal gas accretion rate of the plan- 
ets is calculated. Both models use the criterion that the gas ac- 
cretion rate given by the planet's Kelvin-Helmholtz timescale 
(which we implicitly obtain by the structure calculations) must 
be smaller than the mass transfer rate in the disk. Ida & Lin 
(2004a ) however use this criterion only if the mass of the planet 
is additionally larger than the local gas isolation mass calcu- 
lated from the unperturbed disk profile. This quantity is usually 
clearly larger (Ida & Lin 2004a) than the minimal planet mass 
required to fulfill the viscous and the thermal condition (Lin & 
Papaloizou ll985b . At these masses, the planet becomes able to 
tidally open a (partial) gap in the protoplanetary disk (e.g. Ida 
& Lin 120081) . This could reduce the amount of gas which is ef- 
fectively in the planet's direct gravitational reach, i.e. its gas iso- 
lation mass. Therefore, it seems possible that the planet mass 
where the accretion rate in the disk becomes the limiting fac- 
tor might be significantly smaller than the gas isolation mass as 
assumed in Ida & Lin (2004a). so that we disregard the second 
additional criterion, and always limit the planet's gas accretion 
rate to M disk . 

As illustrated by the prototypical planet of the "main clump", 
the time at which planets start runaway gas accretion occurs for 
typical initial conditions generally at quite advanced stages of 
disk evolution. By then, the disk has undergone significant mass 
loss, and Mdisk is low. Therefore, disk limited planetary accre- 
tion rates are usually down to a few 10~ 4 M e /yr so that growing 
a Jupiter-mass object requires an amount of time comparable to 
the remaining disk lifetime (see § |5.3.2 for a population with very 
long, unrealistic disk lifetimes). Hence, the probability that the 
disk disappears while the planet is at an intermediate mass is not 
vanishingly small. This populates the "planetary desert" with in- 
termediate mass objects, so that we predict only a moderate de- 
crease of the relative number of planets with masses between 30 
to 100 M m (about a factor 2-3 compared to Jovian planets, see the 
planetary initial mass function in paper II). We have investigated 
this point by synthesizing a test population where we limit the 
gas accretion rate of a planet by Mdisk only if its envelope mass 
is larger than the gas isolation mass, as Ida & Lin (2004a), and 
found that this results in a population that has indeed a stronger 
depletion of intermediate mass planets, which are then about 5 
to 6 times less frequent than Jovian planets. 

Long baseline, high precision RV surveys will test for the ex- 
istence of the "horizontal branch" and the "planetary desert", and 
determine their extension and intensity. The upper mass limit of 
the branch constrains the minimal mass needed to go into run- 
away gas accretion, while the lower mass boundary of the branch 
constrains the speed of type I migration. The degree of a deple- 



tion of intermediate mass planets constrains planetary accretion 
rates during runaway gas accretion. 

In the past 12 years, a significant number of extrasolar plan- 
ets with masses much larger than one Jupiter mass were discov- 
ered. We also find such planets (M > 10 Mo) in our model, 
mainly in the "outer group", as shown by fig. 13 Their forma- 
tion is the consequences of the fact that we do not reduce gas 
accretion due to gap formation. About 0.4 % of all initial condi- 
tions lead to the formation of planets with a mass exceeding 12- 
13 M^, which is commonly regarded as the brown dwarf limit 
(Chabrier et al. 120001 . A handful of planets (out of the -50000) 
even reach a mass between 20 to 40 M^. As it was shown re- 
cently that such objects burn deuterium in the layers above the 
core (Baraffe et al. 2008), they form an interesting population of 
deuterium burning planets. 

The synthetic population does not contain analogs of Uranus 
and Neptune both in term of mass and semimajor axis. Indeed, 
no planets at all are found outside ~ 20 AU. Partially, this is sim- 
ply an artifact of the lack of planet formation after the dispersion 
of the gas disk in our model(5 2.4.2 1. Uranus and Neptune are 
however insofar challenging that these planets accreted signifi- 
cant hydrogen-helium envelopes i.e. that they were apparently 
formed while the gas disk was still present (Goldreich et al. 
I2004al Chambers G&Q6]). 

It is well known that core accretion requires for the forma- 
tion of Uranus and Neptune at their current position formation 
timescales which exceed typical gas disk lifetimes by a large 
factor (Pollack et al. 1 19961 Thommes et al. 2003 ). Hence, it was 
proposed that for Uranus and Neptune, one must give up the 
principle of a in-situ formation and take into account the possi- 
bility of an ejection because of N-body interactions with the gi- 
ant planets (Thommes et al. 2002 ; Tsiganis et al. 2005 ). The fact 
that our population contains in the "horizontal branch" a signif- 
icant number of planets with a mass and composition as the ice 
giants, but just at smaller orbital distances, goes well this sec- 
ond interpretation, and the general idea that planetary systems 
start with more crowded and compact configurations (e.g. Ford 
& Chiang[2007]). 



5.3. Non-nominal populations 

The representative tracks we have identified in the previous sec- 
tions and the final a - M were all calculated for one particu- 
lar population, i.e. a particular choice of underlying assumptions 
and parameters of the model like a or the type I migration effi- 
ciency factor. Here we discuss the effect of changing some set- 
tings. 



5.3.1. Core growth regime 

The accretion rate of planetesimals is crucial not only for the 
growth of the cores themselves, but indirectly also for the ac- 
cretion of the envelopes. The solid accretion rate is a function 
of the velocity dispersion Vdi sp of the planetesimals. We use 
the prescription of Pollack et al. (1996]( for Vdi sp corresponding 
to a situation between the shear and the dispersion dominated 
regime. Different results have been obtained concerning the im- 
portance of these two regimes (Ida & Makino I1993I Rafikov 
2003; Goldreich et al. 2004b). In the dispersion dominated pic- 
ture of Thommes et al. (2003 ) planetesimal random velocities 
are higher than in the Pollack et al. (1996) description. Fortier 
et al. (120071) have studied the effects of these high Vd; sp on giant 
planet formation and found an increase of the formation time of 



C. Mordasini et al.: Extrasolar planet population synthesis I 



19 




10* 



1000 



100 



10 




]()♦ 



1000 



100 



10 




1 \i ht liJaiMi ill 



0.1 



a [AU] a [AUJ 

Fig. 14. Planetary formation tracks as in fig.[i]but for non-nominal populations. Top left: For planetesimal eccentricities and incli- 
nations as in Thommes et al. ( 12003 l l. Top right: For M w = (no photoevaporation). Bottom left: For a type I migration efficiency 
factor fi-Q. 1 . Bottom right: For a type II migration rate as in Ida & Lin ( 2004a). 



Jupiter by around one order of magnitude relative to Pollack et 
al. (fl996 ), depending on the mass of the disk. 

We have synthesized a population where we follow Fortier 
et al. ( 120071) in calculating the planetesimal eccentricity and in- 
clination as in Thommes et al. (2003 their eq. 10). The top 
left panel of fig. [14] shows the resulting formation tracks and 
makes clear that there is a very strong effect, especially at large 
distances. Most seeds do not even grow from 0.6 to 1 M®. 
Only an extremely small fraction of initial conditions (extremely 
high fo/c fdisk) allows the formation of giant planets which 
are all inside 0.5 AU and have low masses. It is obvious that 
Kolmogorov-Smirnov tests (paper II) indicate a negligible prob- 
ability that this synthetic population and the observed one come 
from the same parent distribution. We therefore conclude that in 
order to reproduce the observed extrasolar planets by the core 
accretion mechanism, solid accretion must occur on timescales 
clearly shorter than predicted by the model of Thommes et 



al. (|2003). Indeed, including additional effects like planetesi- 
mal fragmentation (Chambers 2006 1 Kenyon & Bromley 2009) 
leads to significantly reduced core growth and planet formation 
timescales (see e.g. Thommes et al. |2008). These effects will 
have, at least qualitatively, similar consequences as the low Vdi sp 
we assume. 



5.3.2. Photoevaporation 

Another population was synthesized with photoevaporation 
switched off (M w = 0), so that the disks evolve only due to 
viscosity. As expected, disk lifetimes are clearly increased for 
such a case. The mean disk lifetime is now about three times 
larger than in the nominal population. This is clearly incompat- 
ible with the observed distribution of disk lifetimes (§|4.3|l. The 
top right panel of fig. 



14 shows that the final positions of the 



planets in the M w = population differ significantly from the 
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nominal case, even if the tracks leading there are similar. As ex- 
pected, migration becomes more important, and a new group of 
massive planets inside about 0.5 AU is formed. The fraction of 
embryos that reach the inner boarder of the computational disk 
increases by about a factor three compared to the nominal case. 
The fraction of initial conditions that leads to giant planets in- 
creases too, as planets have more time to grow which also ex- 
plains why a much emptier "planetary desert" is seen than in 
the nominal case. Statistical tests (paper II) show that this syn- 
thetic population is not compatible with the observed one: The 
KS tests for the Msin/ and the semimajor axis distributions in- 
dicate a probability of only a few percent that the two samples 
come from the same parent population: The synthetic planets 
are too massive and too close-in. We therefore conclude that a 
distribution of disk lifetimes that is incompatible with the obser- 
vational data (Haisch et al. 1200 11 1 leads to a synthetic population 
that is incompatible with the observed exoplanets, too. 



5.3.3. Type I migration 

Motivated by the short migration timescales found by Tanaka et 
al. (120021 ) which indicate that cores are lost to the star before they 
can grow large enough to trigger runaway gas accretion, several 
authors have been re-examining type I migration rates and found 
considerable cause for uncertainties (e.g. Nelson & Papaloizou 
2004, Paardekooper & Mellema 2006). Even if it is clear that 
the simple multiplication of the migration rates of Tanaka et al. 
(2002 ) by a constant factor fi\ is only a first order approximation 
at best of the true type I migration rate (Kley & Crida 2008 ), the 
Monte Carlo simulations represent a suitable way to study the 
effects of various magnitudes of type I migration on the popula- 
tion as a whole, thereby possibly ruling out some values for fi\ 
(Ida & Lin|2008]). 

We have therefore synthesized populations using also 
/i=0.01, 0.1 and 1.0, besides the value of 0.001 used above. In 



14 bottom left, planetary formation track are plotted for a 



population using /j = 0.1 and otherwise nominal settings. For 
the high migration rates, one should keep in mind that migration 
is neglected before the seed embryo is put into the disk thereby 
the true migration is partially underestimated (Ida & Lin 2008 ). 

At fi = 0.1, migration causes many "failed cores" to migrate 
into the feeding limit, especially at small distances as found also 
by Ida & Lin (2008 ). A significant number of "failed cores" nev- 
ertheless remains at intermediate distances, due to the disappear- 
ance of the nebula. Varying the migration rate has also related, 
important consequence on the occurrence of close-in (< 0.3 AU) 
planets with masses 3 < M < 10 M m . Such planets do not ex- 
ist for a low migration rate (cf. figs. [8] and [T3|l but become in- 
creasingly more numerous and smaller with increasing type I 
rate. This absence is the result of the combination of two ef- 
fects. First, our model does not include the in-situ formation of 
terrestrial planets after disk dissipation. Hence, the formation of 
Mercury type planets and larger counterparts, or other mecha- 
nisms that lead to the formation of close-in terrestrial planets 
(see Raymond et al. 120071 for an overview) are not taken into 
account. This in-situ accretion would fill the empty area "from 
below" (in mass). Due to the limited amount of solids very close 
to the star, the mass to which planets can grow in-situ is how- 
ever limited. Second, in most disks seed embryos can only start 
beyond the iceline (cf. the distribution of a sta rt, £ 4.4 1. Thereafter, 



they must be brought in by migration without growing too much 
in order to have a small final mass. This is only possible for fast 
migration (compared to accretion), and sets a minimal mass of 



planets that fill the empty region "from outside" (in semimajor 
axis). 

For fi = 0.1 and 1 .0, the mass of planets that migrate to a to uch 
is usually larger than about 6 and 2 M ffl , respectively. For all fi, 
a few planets with even smaller masses are also found inside 
the feeding limit which come from very small initial positions 
flstart- It should further be noted that this result depends on the 
assumed radial surface density profile of solids. The E D oc cr 3 / 2 
we use for simplicity down to 0. 1 AU leads to high solid surface 
densities near the star. For a truncated disk profile used in tests, 
which might be more appropriate, and fi = 0.1, even planets as 
small as ~ 2 M e have reached the feeding limit, and the region 
inside ~ 0.3 AU which is empty at low fi is then filled. 

Observations able to detect planets in the range 1 to 10 Earth 
masses close to the star will show if there is a minimal mass 
for planets inside a few tens of an AU, and whether the mass 
spectrum there is continuous, or contains a gap between plan- 
ets formed in-situ and cores brought in from further out by mi- 
gration. Hence, such object will potentially provide us with an 
indicator on the efficiency of type I migration. 

Considering the very massive planets (> 10 M^), it is found 
that at low type I migration rates, no such planets are found 
within about 0.4 AU of the star. At high rates, some very massive 
planets have, in contrast, migrated to the feeding limit. It is in- 
teresting to note that the type I migration, which affects directly 
small planets, also has a significant and measurable effect on the 
most massive planets as it allows planets forming at small dis- 
tances to have a quicker access to larger amounts of solids than 
the locally available M iso for fi = 0.001. 

5.3.4. Type II migration 

The formation tracks of giant planets of both of the "main 
clump" and the "outer group" show the importance of the planet- 
dominated migration phase, as well as the disk limited gas accre- 
tion rate (which itself results from the process of gas transport 
inside the protoplanetary disk) in the braking phase. It is inter- 
esting to note that, using planet formation models following the 
growth and migration of a large number of embryos, Thommes 
et al. d20081 > came to the same conclusion regarding the impor- 
tance of these two processes. 

The braking phase plays a crucial role for the final a-M posi- 
tion of gas rich planets. It is clear that our description of the plan- 
etary dm/dt and da/dt in this phase remains a first approxima- 
tion that could be significantly improved by a more sophisticated 
description of planet-disk interactions. For example, the eccen- 
tric instability that justifies the assumption that the planetary gas 
accretion rate is the same as the disk accretion rate occurs only 
for planets larger than a certain mass (Kley & Dirksen 2006|. At 
smaller masses, the accretion rate could be smaller (D' Angelo et 
a U20021 l. Also planet-planet interactions of several giant planets 
forming concurrently can have important consequences on the 
migration and accretion rate. This was shown by Thommes et 
al. ©OS), where a chain of migrating giant planets, locked to- 
gether by mean motion resonances, is a frequent configuration. 
Then, only the outermost planet directly exchanges torques with 
the gas disk, and accretes gas from it. While with our one embryo 
per disk approach, we obviously cannot model such a behavior, 
we note that this configuration typically leads to an even earlier 
transition to the planet-dominated regime as multiple planets are 
better at holding off the outer disk than one (Thommes et al. 
120081) . 

The da/dt of a satellite that is massive compared to the local 
disk mass was studied by Syer & Clarke (1995) and Ivanov et 
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al. ©99). They both agree that in this situation, the migration 
timescale becomes longer than the viscous timescale, but quanti- 
tatively their results differ. The differences are linked to the ques- 
tion of what part of the angular momentum flux is effective in 
moving the planet (Ida & Lin 2005|, and whether gas overflows 
the gap formed by the planet (Armitage 2007). Therefore, there 
is some uncertainty on the exact migration rate in the case of a 
massive planet. 

To have an estimate of the sensitivity of our results on the mi- 
gration rate in the braking phase we have synthesized a test pop- 
ulation replacing eq.[3]by the prescription of Ida & Lin (2004a ). 



The bottom right panel of fig. 14 shows the resulting forma- 
tion tracks. One sees that there is a visible modification of the 
tracks, but that the basic types remain. The comparison of the 
tracks with the dashed lines (slope equal —n) shows that with 
this alternative prescription, type II migration occurs on a shorter 
timescale. The reason is that the da/dt is now proportional to 
I.(R m )Rl instead of I,(a pi . dnet )a^ liineV where R,„ is the radius of 
maximum viscous couple, which moves outwards very quickly, 
and is in almost all cases larger than a p i an et once planets are large 
enough to migrate in type II. Therefore is the first quantity usu- 
ally larger than the second one (fig. [2]). 

There are also two examples of tracks where outward migra- 
tion occurs. As expected from the quick outward movement of 
R m , only a very small group of about 10 seeds (out of ~ 10000 
initial conditions) has undergone significant outward migration^] 
The maximal planetary mass M max at ~ 20 AU is therefore in- 



creased relative to the nominal model (sect. 5.2.2 1 by about one 
order of magnitude. 

The discovery of many more giant planets outside ~ 20 AU 
by techniques like astrometry or direct imaging (Kalas et al. 
2008 ; Marois et al. 2008 ) would indicate that besides outward 
migration one or several of the following mechanisms not in- 
cluded in the nominal model is important: Scattering of plan- 
etary seeds to large a start early during formation ("monarchical 
growth", Weidenschilling 2005, 2008 ), a completely different 
formation mechanism (direct gravitational collapse, e.g. Boss 
120011 Mayer et al. 120041 ) or scattering of planets in initially 
crowded multiple planetary systems after formation (e.g. Veras 
& Armitage 2 003t . Note however for the last possibility that the 
simulations of Thommes et al. (12008 ), which combine in a self- 
consistent way the formation and subsequent N-body interac- 
tions, do not usually lead to the formation of giant planets further 
out than ~ 20 AU. In another non-nominal population, we have 
addressed the early embryo ejection and used a modified pre- 
scription to calculate the starting time f start that approximately 
mimics a "monarchical growth" mode (Weidenschilling 2005 
2008 ). In this case, while leaving the population in the inner sys- 
tem relatively untouched, a population of very massive planets 
outside 3 AU comes into existence so that the maximal plane- 
tary mass M max does not decrease anymore outside ~ 3 AU, but 
remains constant out to ~ 10 - 20 AU. 



5.3.5. Transition from type I to type II migration 

As mentioned at the beginning of the paper, the transition from 
type I to type II occurs when the planet can open a gap in the 
protoplanetary disk. In our nominal model we assume that this 
occurs when the disk scale height becomes smaller than the 



3 Note that we calculate in this test population R,„ as Ida & Lin 
(2004a), and not in a self-consistent way from the disk model, where 



photoevaporation can influence the evolution of R,„ at late times and 
hence the outward migration (Veras & Armitage 2003 1. 



planet's Hill radius (which is the relevant criterion in the limit 
of very large Reynolds numbers). In order to infer the influence 
of this part of our model, we have synthesized non-nominal pop- 
ulations using the type I/type II transition condition derived by 
Crida et al. (2006). Using the low type I migration of the nominal 
model, only very few "Hot" planets are then formed, since the 
transition to (faster) type II occurs at higher masses. However, 
when fi is increased, we obtain similar formation tracks and fi- 
nal sub-populations. Note that in this case, very few planet mi- 
grate in the disk-dominated type II regime, and a majority of 
them passes directly from type I to the planet-dominated type 
II regime. We mention finally that, as outlined by e.g. Armitage 
and Rice (2005 ), the transition, and especially the migration rate 
during the transition, is not fully understood and could result in 
special migration regimes not considered here (e.g. Papaloizou 
& Terquem 2006). This incomplete understanding of migration 
in general remains one of the major uncertainties in planet for- 
mation theories of today. 

6. Summary and conclusion 

We have presented our extrasolar planet population synthe- 
sis calculations. As formation model we use a slightly simpli- 
fied version of the extended core accretion model presented in 
Alibert et al. j2005al ) which has been shown to reproduce many 
observed properties of the giant planets of our own solar system 
(Alibert et al.|2005bj. 

We use four random variables to describe the possible initial 
conditions for planet formation: The dust-to-gas ratio /d/Gi the 
initial gas surface density £o, the photoevaporation rate M w and 
the starting position of the embryo in the disk, a star t. The distri- 
butions for the first three Monte Carlo variables can be derived 
from observed properties of stars of the solar neighborhood, or 
protoplanetary disks. For the dust-to-gas ratio, we use the dis- 
tribution of [Fe/H] of the stars in the CORALIE planet search 
sample (Santos et. 2003). For the gas surface densities, we use 
the disk mass distribution in p Ophiuchi (Beckwith & Sargent 
1996). For the photoevaporation rate, we use the distribution of 
disk lifetimes by Haisch et al. (2001 ). We have also discussed the 
complications that arise from the conversion of such observed 
quantities into figures that can be used in numerical simulations. 
The last random variable, a start cannot be derived from observa- 
tions. Here we follow Ida & Lin (2004a ) and use a distribution 
that is uniform in log(a). 

With the adopted random distributions, we find that our for- 
mation model predicts a population of synthetic planets that is 
very diverse, not unlike the actually observed population. The 
final semimajor axis varies by two orders of magnitude, and the 
mass by four orders. Similarly, the internal bulk composition is 
very different, and covers gas giants but also "Super Earth" mass 
planets with a envelope to core mass ratio similar to Venus. We 
conclude that this diversity illustrates the ability of the underly- 
ing concepts of core formation to serve as an unified formation 
model for planets of masses ranging from a few times the Earth 
mass to beyond the brown dwarf limit. The observed diversity of 
extrasolar planet is a natural consequence of the different prop- 
erties of protoplanetary disks. 

Planet formation is a process of concurrent mass accretion 
and migration which can be well represented by planetary for- 
mation tracks in the mass-distance plane. These tracks show a 
number of distinct phases, which are visible in the final a - M 
distribution of the planets, too. 

In the first phase, at small masses, planets migrate in type I 
migration (which must however be slow to be compatible with 
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observations, cf. paper II), and accrete mostly solids. Planets 
which remain in this stage until the moment when the disk dis- 
appears form a vast sub-population of low mass, core-dominated 
planets, "failed cores". The distributions of protoplanetary disks 
properties are such that this is the most likely outcome of the 
formation process. This is consistent with the observation that 
90 to 95% of FGK stars in the solar neighborhood apparently 
remained without giant planets. 

If the disk properties are instead such that the planet grows 
massive enough to open a gap in the disk, the second phase 
starts. Then, planets migrate inwards in type II migration, col- 
lecting solids on their way in. In this phase, migration occurs on 
a shorter timescale than accretion, so that the tracks of the plan- 
ets show a "horizontal branch" phase in the a - M plot. Planets 
in this phase have masses between 10 to 30 M s , a static gaseous 
envelope and an internal gas to solid ratio similar to Neptune. 

Some planets on the "horizontal branch" collect enough 
solids to become supercritical for gas runaway accretion. Their 
gas accretion rate then quickly reaches the disk limited value, 
and their mass becomes large compared to the local disk mass. 
Therefore, giant gaseous planets reach their final mass and semi- 
major axis in the third phase, when they accrete gas at the rate 
controlled by the disk and migrate in the planet dominated, 
slower type II regime. This leads to the formation of a concentra- 
tion of giant planets between 0.3 to 2 AU in the "main clump", 
with masses between ~ 100 M e and 3 Mq, . 

Embryos starting far from the parent star in metal rich disk 
can grow supercritical for gas runaway accretion in-situ, with- 
out passing through the "horizontal branch". This leads to the 
formation of an "outer group" of giant planets. 

Apart from these four sub-populations, we find in the final 
nominal a - M diagram (1) an absence of massive (> lOM^) 
planets both close (< 0.5 AU) to the star and very far (> 10 AU) 
from it, (2) a certain depletion of planets with masses between 
30 and 100 M m , in analogy to the "planetary desert" (Ida & Lin 
2004a) which is however not very strong (a factor 2-3 relative to 
giant planets), as at the time planets start runaway gas accretion, 
the disk limited accretion rate, which is decisive in this regime, is 
usually already quite low, (3) a handful very massive (> 20Mt ) _) 
deuterium burning planets (Baraffe et al. 2008), showing that this 
mass domain can be populated by different formation channels, 
(4) an absence of Neptune analogs in terms of both mass and 
semimajor axis, but many such planets at smaller distances. 

Population synthesis is a valuable tool to asses the global 
consequences of model settings. From the synthesis of a num- 
ber of non-nominal populations we find that (1) cores must form 
quickly (Thommes et al. 120031 Chamber |2006), (2) disks with 
lifetimes too long compared to observations lead to too massive 
planets too close-in, (3) fast type I migration brings low mass 
planets close to the parent star, so that the depletion of planets 
with 3 < M/M e < 10 inside 0.3 AU in the nominal popula- 
tion vanishes, (4) the early ejection of seeds, and outward type 
II migration could result in very massive planets out to about 20 
AU, but only in small numbers, (5) the transition regime between 
type I and type II migration has important implications for the 
final population. 

With improving detection methods leading to a rapidly in- 
creasing number of known extrasolar planets, the actual distri- 
butions of masses and semi-major axis can be determined with 
growing statistical confidence and compared to model predic- 
tions. Hence, planet populations studies are the ideal tools to 
extract from these distributions constraints on various aspects of 
the formation models thereby reaping the scientific benefits of 
the large observational efforts invested. In the companion paper 



(paper II), we carry out a direct statistical comparison between 
a carefully chosen sample of actual exoplanets and detectable 
synthetic planets and use this comparison to extract relevant con- 
straints on the formation mechanism. Last but not least, popula- 
tion synthesis can itself be used to guide and perhaps optimize 
the next generation of detection instruments. 
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